1 Exploring the data set

The HS data set was previously used in the CEI paper (Martenies et al., 2019). In the original analysis, we used an exposure index based on the CalEnvironScreen tool. We observed lower birth weights and lower adiposity associated with higher index scores, driven largely by exposures to social indicators of health at the neighborhood level. Now, we are aiming to use methods for mixtures to try to identify which exposures are driving these association.

The complete data set for the birth weight outcome consists of n = 897 participants. This represents 63.62% of the original Healthy Start 1 cohort.

Of the 897 participants, 0.27% identify as Latina, 0.17% identify as Black, and 0.27% identify as another non-NHW race or ethnicity. The median age of mothers in this dataset is 28 years. 0.51% of babies born were male.

1.1 Exposure data

We have included 19 exposures in our analysis.

These exposures are based on the census tract where each mother lived at the time of enrollment into Healthy Start. With the exception of air pollution (mean_pm and mean_o3), these are based on long-term averages at for each census tract. For mean_pm and mean_o3 are based on the average pollution levels across each pregnancy (est. conception date to delivery date) estimated using ordinary kriging and monitoring data.

#' Exposure data
X <- select(hs_data2, mean_pm, mean_o3, pct_tree_cover, pct_impervious,
            mean_aadt_intensity, dist_m_tri:dist_m_mine_well,
            cvd_rate_adj, res_rate_adj, violent_crime_rate, property_crime_rate,
            pct_less_hs, pct_unemp, pct_limited_eng, pct_hh_pov) %>%
  as.matrix()
head(X)
##       mean_pm  mean_o3 pct_tree_cover pct_impervious mean_aadt_intensity
## [1,] 8.483046 47.19072       6.006276       43.30893          10128.4962
## [2,] 6.598608 50.05090       7.281109       48.36432          10749.0359
## [3,] 7.454146 48.57052      17.205991       31.67281           9048.6468
## [4,] 6.671239 50.06429       6.842898       45.00359           4223.3434
## [5,] 7.122537 50.14275       3.357792       28.16745            858.7283
## [6,] 7.637453 47.03125      10.743612       45.87564          15603.9800
##      dist_m_tri dist_m_npl dist_m_waste_site dist_m_major_emit dist_m_cafo
## [1,]   2827.538   729.2371          4829.780          7968.654    29116.58
## [2,]   1576.420  5239.2211          4417.792          3780.951    51044.30
## [3,]   3350.303  2992.2968          5211.871          7423.232    36079.21
## [4,]   3364.954  6998.1286          8921.318          9636.816    42235.78
## [5,]   2923.811  3427.2247          7006.042          6806.912    29145.98
## [6,]   3364.200  3166.5395          4484.960          5265.285    43921.85
##      dist_m_mine_well cvd_rate_adj res_rate_adj violent_crime_rate
## [1,]        1749.1256     275.2480     155.7767          14.377133
## [2,]        7354.5310     279.6435     226.8038           8.905404
## [3,]        4887.2996     221.0414     157.6974           7.636888
## [4,]        3752.6399     203.8812     142.5368           2.850212
## [5,]         729.7784     194.1983     101.0046           5.435988
## [6,]        5870.6867     174.3361     120.3281           5.035971
##      property_crime_rate pct_less_hs pct_unemp pct_limited_eng pct_hh_pov
## [1,]            37.32935   31.784946 11.529628       26.114650  12.010919
## [2,]            67.03932   15.290231  4.908306        8.500401  18.123496
## [3,]            46.78194    6.891702  4.564963        0.000000   6.307978
## [4,]            21.95270    2.725915  5.623583        1.350621   9.292274
## [5,]            22.49834   12.919186  5.234103        6.307385   2.115768
## [6,]            47.15500    3.842365 10.000000        5.121799  25.171768

Variance and histograms of the exposure variables (in their original units):

var(X)
##                             mean_pm        mean_o3 pct_tree_cover
## mean_pm                 0.391784015    0.006083605     -0.2054297
## mean_o3                 0.006083605    9.383489039     -0.4158151
## pct_tree_cover         -0.205429726   -0.415815089      9.7193077
## pct_impervious          0.508898445   -1.674151031      5.8719893
## mean_aadt_intensity  -182.234953786  474.627052967   8431.6446632
## dist_m_tri           -255.176839682  444.286548683    -73.1423054
## dist_m_npl           -289.002141382  539.849185829    434.4654007
## dist_m_waste_site    -275.262105884  261.902915064   1933.8647304
## dist_m_major_emit      71.096593638  577.257325397    265.4284518
## dist_m_cafo         -1291.237441927  -35.275020052  10170.6234275
## dist_m_mine_well     -339.250592215 -375.434990683   3136.3680766
## cvd_rate_adj            3.871688575    0.939328342    -24.8232924
## res_rate_adj            2.356328835   -0.181515705     -3.5331376
## violent_crime_rate      0.232839920    0.577648302     -4.0583754
## property_crime_rate     2.001989749   -2.773092354    -22.6429724
## pct_less_hs             1.132232861    0.637361326     -7.5753471
## pct_unemp               0.100439902    0.288530482     -0.3330523
## pct_limited_eng         0.432516169    0.295617023     -2.8349116
## pct_hh_pov              0.731824476   -0.606648513      0.3805472
##                     pct_impervious mean_aadt_intensity     dist_m_tri
## mean_pm                  0.5088984           -182.2350     -255.17684
## mean_o3                 -1.6741510            474.6271      444.28655
## pct_tree_cover           5.8719893           8431.6447      -73.14231
## pct_impervious         176.8316214          55459.6063   -15279.44024
## mean_aadt_intensity  55459.6063235       67283287.0201 -1315386.69307
## dist_m_tri          -15279.4402428       -1315386.6931  6558190.20296
## dist_m_npl           -7729.3843793        1683196.0799  4282727.94125
## dist_m_waste_site    -4662.9983638        2039577.9230  2441267.84540
## dist_m_major_emit     2627.0270993        2477155.3406  1433153.16531
## dist_m_cafo          16586.9964129       15462371.9832  3431065.70215
## dist_m_mine_well       706.6674650        2073244.5987   995872.11873
## cvd_rate_adj           230.4542985          20477.4374   -49347.60273
## res_rate_adj           176.8108084          33055.3733   -31870.98664
## violent_crime_rate      26.6945028           5736.5627    -1014.08753
## property_crime_rate    118.0737725          22077.3894    -5365.69285
## pct_less_hs             56.8383947          -4056.6889   -12372.14262
## pct_unemp               25.9434246           6003.3343    -2527.22451
## pct_limited_eng         41.9919053           2620.6198    -5408.86434
## pct_hh_pov              82.2198624          17850.1649    -8842.76408
##                        dist_m_npl dist_m_waste_site dist_m_major_emit
## mean_pm                 -289.0021         -275.2621          71.09659
## mean_o3                  539.8492          261.9029         577.25733
## pct_tree_cover           434.4654         1933.8647         265.42845
## pct_impervious         -7729.3844        -4662.9984        2627.02710
## mean_aadt_intensity  1683196.0799      2039577.9230     2477155.34057
## dist_m_tri           4282727.9413      2441267.8454     1433153.16531
## dist_m_npl          11125411.7474      4193498.0586     6948817.25739
## dist_m_waste_site    4193498.0586      5344101.7540     1395277.06805
## dist_m_major_emit    6948817.2574      1395277.0681    10114549.72263
## dist_m_cafo          5416531.1320      5586018.8251    -2993791.05377
## dist_m_mine_well      256924.3029      1375784.7856    -1810174.74785
## cvd_rate_adj          -30921.0390       -43119.5785       16272.40152
## res_rate_adj          -19393.1304       -32402.8440       -1320.21297
## violent_crime_rate      -672.9264        -3702.6112        -360.49700
## property_crime_rate   -18283.4264       -22350.3006      -24007.42305
## pct_less_hs            -6760.5337       -11422.4985        8866.74917
## pct_unemp               2195.0515        -1476.4094        5212.74830
## pct_limited_eng          498.0033        -4277.8134        9367.28435
## pct_hh_pov             -1135.3843        -7599.7432        8682.26135
##                        dist_m_cafo dist_m_mine_well   cvd_rate_adj
## mean_pm                -1291.23744        -339.2506      3.8716886
## mean_o3                  -35.27502        -375.4350      0.9393283
## pct_tree_cover         10170.62343        3136.3681    -24.8232924
## pct_impervious         16586.99641         706.6675    230.4542985
## mean_aadt_intensity 15462371.98316     2073244.5987  20477.4373759
## dist_m_tri           3431065.70215      995872.1187 -49347.6027339
## dist_m_npl           5416531.13199      256924.3029 -30921.0389720
## dist_m_waste_site    5586018.82514     1375784.7856 -43119.5785165
## dist_m_major_emit   -2993791.05377    -1810174.7478  16272.4015197
## dist_m_cafo         46324000.89481     9345575.3772 -46645.9665229
## dist_m_mine_well     9345575.37722     4430024.9964 -39046.5984701
## cvd_rate_adj          -46645.96652      -39046.5985   2039.8569530
## res_rate_adj          -13772.40263      -16322.5110   1289.5661935
## violent_crime_rate       722.31907       -2032.3464    135.9487143
## property_crime_rate   -15833.92381       -4272.3829    343.9364726
## pct_less_hs           -26060.83378      -10037.6577    328.3044447
## pct_unemp              -1030.96916       -2827.2369    105.0153846
## pct_limited_eng        -7089.15821       -4814.6687    183.5853966
## pct_hh_pov              -855.38016       -5030.4055    266.1004715
##                       res_rate_adj violent_crime_rate property_crime_rate
## mean_pm                  2.3563288          0.2328399            2.001990
## mean_o3                 -0.1815157          0.5776483           -2.773092
## pct_tree_cover          -3.5331376         -4.0583754          -22.642972
## pct_impervious         176.8108084         26.6945028          118.073773
## mean_aadt_intensity  33055.3733277       5736.5627383        22077.389365
## dist_m_tri          -31870.9866403      -1014.0875345        -5365.692846
## dist_m_npl          -19393.1304345       -672.9263612       -18283.426420
## dist_m_waste_site   -32402.8439544      -3702.6111771       -22350.300554
## dist_m_major_emit    -1320.2129699       -360.4970006       -24007.423046
## dist_m_cafo         -13772.4026269        722.3190727       -15833.923813
## dist_m_mine_well    -16322.5110008      -2032.3464340        -4272.382880
## cvd_rate_adj          1289.5661935        135.9487143          343.936473
## res_rate_adj          1091.1856742        104.4979610          333.780710
## violent_crime_rate     104.4979610         40.1175363          160.725724
## property_crime_rate    333.7807097        160.7257236         1295.004010
## pct_less_hs            197.8827546         22.5579950           -3.138375
## pct_unemp               72.3576933         11.3130282            1.362247
## pct_limited_eng        104.0524036         12.7978322          -14.963510
## pct_hh_pov             201.6582659         29.1947400           64.236239
##                        pct_less_hs     pct_unemp pct_limited_eng    pct_hh_pov
## mean_pm                  1.1322329     0.1004399       0.4325162     0.7318245
## mean_o3                  0.6373613     0.2885305       0.2956170    -0.6066485
## pct_tree_cover          -7.5753471    -0.3330523      -2.8349116     0.3805472
## pct_impervious          56.8383947    25.9434246      41.9919053    82.2198624
## mean_aadt_intensity  -4056.6889048  6003.3343312    2620.6197529 17850.1649192
## dist_m_tri          -12372.1426191 -2527.2245090   -5408.8643368 -8842.7640785
## dist_m_npl           -6760.5337115  2195.0514738     498.0033420 -1135.3843390
## dist_m_waste_site   -11422.4985495 -1476.4094188   -4277.8133935 -7599.7432386
## dist_m_major_emit     8866.7491706  5212.7483023    9367.2843472  8682.2613524
## dist_m_cafo         -26060.8337755 -1030.9691591   -7089.1582114  -855.3801591
## dist_m_mine_well    -10037.6576614 -2827.2368665   -4814.6687400 -5030.4055237
## cvd_rate_adj           328.3044447   105.0153846     183.5853966   266.1004715
## res_rate_adj           197.8827546    72.3576933     104.0524036   201.6582659
## violent_crime_rate      22.5579950    11.3130282      12.7978322    29.1947400
## property_crime_rate     -3.1383751     1.3622468     -14.9635105    64.2362387
## pct_less_hs            162.1681017    39.4206217      85.1910014   100.9072175
## pct_unemp               39.4206217    24.6546969      25.2172769    36.9693212
## pct_limited_eng         85.1910014    25.2172769      68.6532943    67.2758215
## pct_hh_pov             100.9072175    36.9693212      67.2758215   119.7157808
ggplot(pivot_longer(as.data.frame(X), mean_pm:pct_hh_pov, names_to = "exp", values_to = "value")) + 
    geom_histogram(aes(x = value)) + 
    facet_wrap(~ exp, scales = "free")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Scaling the exposure variables

X.scaled <- apply(X, 2, scale)
head(X.scaled)
##          mean_pm    mean_o3 pct_tree_cover pct_impervious mean_aadt_intensity
## [1,]  1.60876944 -0.2502907    -0.08261827      0.2084897         -0.02626143
## [2,] -1.40186806  0.6834152     0.32629926      0.5886571          0.04938980
## [3,] -0.03503482  0.2001460     3.50981981     -0.6665506         -0.15790801
## [4,] -1.28583023  0.6877893     0.18573791      0.3359288         -0.74617032
## [5,] -0.56482237  0.7133998    -0.93215008     -0.9301550         -1.15635722
## [6,]  0.25782234 -0.3023500     1.43693690      0.4015075          0.64126566
##      dist_m_tri  dist_m_npl dist_m_waste_site dist_m_major_emit dist_m_cafo
## [1,] -0.4008664 -1.44212141      -0.172775897        -0.1090638  -1.1354079
## [2,] -0.8894134 -0.08999604      -0.350992130        -1.4258114   2.0863310
## [3,] -0.1967327 -0.76363997      -0.007492447        -0.2805619  -0.1124210
## [4,] -0.1910117  0.43733697       1.597126043         0.4154599   0.7921355
## [5,] -0.3632728 -0.63324550       0.768623456        -0.4743524  -1.1310891
## [6,] -0.1913062 -0.71140077      -0.321936759        -0.9590895   1.0398621
##      dist_m_mine_well cvd_rate_adj res_rate_adj violent_crime_rate
## [1,]       -0.7748959    0.6917151   -0.2847192          0.2444833
## [2,]        1.8883051    0.7890362    1.8654606         -0.6194047
## [3,]        0.7160914   -0.5084805   -0.2265752         -0.8196807
## [4,]        0.1769998   -0.8884275   -0.6855261         -1.5754111
## [5,]       -1.2592010   -1.1028189   -1.9428184         -1.1671633
## [6,]        1.1833114   -1.5425892   -1.3578439         -1.2303189
##      property_crime_rate pct_less_hs   pct_unemp pct_limited_eng pct_hh_pov
## [1,]          -0.5008789  1.20091488  0.36793456      2.15383825 -0.3020187
## [2,]           0.3247153 -0.09436047 -0.96557105      0.02798425  0.2566428
## [3,]          -0.2382061 -0.75386917 -1.03471894     -0.99792447 -0.8232412
## [4,]          -0.9281724 -1.08099463 -0.82151747     -0.83491873 -0.5504902
## [5,]          -0.9130100 -0.28055077 -0.89995692     -0.23668961 -1.2063898
## [6,]          -0.2278392 -0.99332353  0.05987409     -0.37977738  0.9008223

Variance and histograms of the exposure variables (scaled):

var(X.scaled)
##                          mean_pm      mean_o3 pct_tree_cover pct_impervious
## mean_pm              1.000000000  0.003172893   -0.105274276     0.06114033
## mean_o3              0.003172893  1.000000000   -0.043541201    -0.04109912
## pct_tree_cover      -0.105274276 -0.043541201    1.000000000     0.14164056
## pct_impervious       0.061140333 -0.041099117    0.141640558     1.00000000
## mean_aadt_intensity -0.035493982  0.018889337    0.329716765     0.50844411
## dist_m_tri          -0.159193706  0.056635532   -0.009161339    -0.44867872
## dist_m_npl          -0.138426634  0.052836279    0.041781063    -0.17426369
## dist_m_waste_site   -0.190232937  0.036984589    0.268331135    -0.15168685
## dist_m_major_emit    0.035715124  0.059253499    0.026770503     0.06211712
## dist_m_cafo         -0.303095664 -0.001691929    0.479321466     0.18326720
## dist_m_mine_well    -0.257510042 -0.058230361    0.477976199     0.02524829
## cvd_rate_adj         0.136954783  0.006789463   -0.176295758     0.38371166
## res_rate_adj         0.113962827 -0.001793836   -0.034307854     0.40251263
## violent_crime_rate   0.058730941  0.029772424   -0.205526308     0.31693831
## property_crime_rate  0.088879773 -0.025156291   -0.201827442     0.24673907
## pct_less_hs          0.142046220  0.016338825   -0.190810449     0.33564418
## pct_unemp            0.032317153  0.018969666   -0.021515177     0.39291400
## pct_limited_eng      0.083396591  0.011647067   -0.109746623     0.38111402
## pct_hh_pov           0.106858205 -0.018100029    0.011156172     0.56509450
##                     mean_aadt_intensity   dist_m_tri  dist_m_npl
## mean_pm                     -0.03549398 -0.159193706 -0.13842663
## mean_o3                      0.01888934  0.056635532  0.05283628
## pct_tree_cover               0.32971677 -0.009161339  0.04178106
## pct_impervious               0.50844411 -0.448678724 -0.17426369
## mean_aadt_intensity          1.00000000 -0.062619247  0.06152095
## dist_m_tri                  -0.06261925  1.000000000  0.50138396
## dist_m_npl                   0.06152095  0.501383960  1.00000000
## dist_m_waste_site            0.10755964  0.412369055  0.54385239
## dist_m_major_emit            0.09495686  0.175965418  0.65505772
## dist_m_cafo                  0.27696155  0.196849357  0.23859436
## dist_m_mine_well             0.12008641  0.184760224  0.03659688
## cvd_rate_adj                 0.05527416 -0.426652406 -0.20525615
## res_rate_adj                 0.12199419 -0.376750794 -0.17601130
## violent_crime_rate           0.11041575 -0.062519618 -0.03185242
## property_crime_rate          0.07479259 -0.058223492 -0.15232247
## pct_less_hs                 -0.03883608 -0.379376300 -0.15916230
## pct_unemp                    0.14739715 -0.198747643  0.13253691
## pct_limited_eng              0.03855846 -0.254907958  0.01801953
## pct_hh_pov                   0.19888999 -0.315587892 -0.03111065
##                     dist_m_waste_site dist_m_major_emit  dist_m_cafo
## mean_pm                   -0.19023294        0.03571512 -0.303095664
## mean_o3                    0.03698459        0.05925350 -0.001691929
## pct_tree_cover             0.26833114        0.02677050  0.479321466
## pct_impervious            -0.15168685        0.06211712  0.183267205
## mean_aadt_intensity        0.10755964        0.09495686  0.276961553
## dist_m_tri                 0.41236906        0.17596542  0.196849357
## dist_m_npl                 0.54385239        0.65505772  0.238594356
## dist_m_waste_site          1.00000000        0.18977973  0.355027509
## dist_m_major_emit          0.18977973        1.00000000 -0.138307324
## dist_m_cafo                0.35502751       -0.13830732  1.000000000
## dist_m_mine_well           0.28275484       -0.27042332  0.652378929
## cvd_rate_adj              -0.41298787        0.11328659 -0.151743889
## res_rate_adj              -0.42432287       -0.01256670 -0.061257230
## violent_crime_rate        -0.25287368       -0.01789622  0.016755559
## property_crime_rate       -0.26866460       -0.20976678 -0.064647236
## pct_less_hs               -0.38800832        0.21893159 -0.300678627
## pct_unemp                 -0.12862329        0.33009857 -0.030506530
## pct_limited_eng           -0.22333346        0.35547555 -0.125707430
## pct_hh_pov                -0.30045944        0.24950766 -0.011486308
##                     dist_m_mine_well cvd_rate_adj res_rate_adj
## mean_pm                  -0.25751004  0.136954783  0.113962827
## mean_o3                  -0.05823036  0.006789463 -0.001793836
## pct_tree_cover            0.47797620 -0.176295758 -0.034307854
## pct_impervious            0.02524829  0.383711656  0.402512635
## mean_aadt_intensity       0.12008641  0.055274160  0.121994187
## dist_m_tri                0.18476022 -0.426652406 -0.376750794
## dist_m_npl                0.03659688 -0.205256148 -0.176011297
## dist_m_waste_site         0.28275484 -0.412987865 -0.424322872
## dist_m_major_emit        -0.27042332  0.113286594 -0.012566704
## dist_m_cafo               0.65237893 -0.151743889 -0.061257230
## dist_m_mine_well          1.00000000 -0.410752544 -0.234765650
## cvd_rate_adj             -0.41075254  1.000000000  0.864359590
## res_rate_adj             -0.23476565  0.864359590  1.000000000
## violent_crime_rate       -0.15245003  0.475234675  0.499449246
## property_crime_rate      -0.05640681  0.211613232  0.280786581
## pct_less_hs              -0.37449548  0.570813439  0.470409304
## pct_unemp                -0.27052616  0.468277441  0.441149256
## pct_limited_eng          -0.27607853  0.490577454  0.380164971
## pct_hh_pov               -0.21843600  0.538480631  0.557944498
##                     violent_crime_rate property_crime_rate pct_less_hs
## mean_pm                     0.05873094         0.088879773  0.14204622
## mean_o3                     0.02977242        -0.025156291  0.01633882
## pct_tree_cover             -0.20552631        -0.201827442 -0.19081045
## pct_impervious              0.31693831         0.246739067  0.33564418
## mean_aadt_intensity         0.11041575         0.074792588 -0.03883608
## dist_m_tri                 -0.06251962        -0.058223492 -0.37937630
## dist_m_npl                 -0.03185242        -0.152322474 -0.15916230
## dist_m_waste_site          -0.25287368        -0.268664603 -0.38800832
## dist_m_major_emit          -0.01789622        -0.209766780  0.21893159
## dist_m_cafo                 0.01675556        -0.064647236 -0.30067863
## dist_m_mine_well           -0.15245003        -0.056406808 -0.37449548
## cvd_rate_adj                0.47523468         0.211613232  0.57081344
## res_rate_adj                0.49944925         0.280786581  0.47040930
## violent_crime_rate          1.00000000         0.705151942  0.27967307
## property_crime_rate         0.70515194         1.000000000 -0.00684836
## pct_less_hs                 0.27967307        -0.006848360  1.00000000
## pct_unemp                   0.35971778         0.007623781  0.62343462
## pct_limited_eng             0.24385889        -0.050184228  0.80738433
## pct_hh_pov                  0.42127121         0.163143151  0.72420883
##                        pct_unemp pct_limited_eng  pct_hh_pov
## mean_pm              0.032317153      0.08339659  0.10685820
## mean_o3              0.018969666      0.01164707 -0.01810003
## pct_tree_cover      -0.021515177     -0.10974662  0.01115617
## pct_impervious       0.392914001      0.38111402  0.56509450
## mean_aadt_intensity  0.147397153      0.03855846  0.19888999
## dist_m_tri          -0.198747643     -0.25490796 -0.31558789
## dist_m_npl           0.132536906      0.01801953 -0.03111065
## dist_m_waste_site   -0.128623290     -0.22333346 -0.30045944
## dist_m_major_emit    0.330098571      0.35547555  0.24950766
## dist_m_cafo         -0.030506530     -0.12570743 -0.01148631
## dist_m_mine_well    -0.270526163     -0.27607853 -0.21843600
## cvd_rate_adj         0.468277441      0.49057745  0.53848063
## res_rate_adj         0.441149256      0.38016497  0.55794450
## violent_crime_rate   0.359717785      0.24385889  0.42127121
## property_crime_rate  0.007623781     -0.05018423  0.16314315
## pct_less_hs          0.623434625      0.80738433  0.72420883
## pct_unemp            1.000000000      0.61293958  0.68048090
## pct_limited_eng      0.612939575      1.00000000  0.74208323
## pct_hh_pov           0.680480902      0.74208323  1.00000000
ggplot(pivot_longer(as.data.frame(X.scaled), mean_pm:pct_hh_pov, 
                    names_to = "exp", values_to = "value")) + 
    geom_histogram(aes(x = value)) + 
    facet_wrap(~ exp, scales = "free")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.2 Covariate data

Covariates were assessed at the individual level. These were selected based on previous HS studies and others in the literature and informed by a DAG.

NOTE: It’ll be interesting to see what comes out of our BEAMERS discussion re: adjusting for gestational age. It’s currently in the analysis

There are four continuous covariates; all of the others have been coded as dummy variables. For the dummy variables, the reference groups are: white_re, ed_grad, norm_bmi

W <- select(hs_data2, latina_re, black_re, other_re,
            ed_no_hs, ed_hs, ed_aa, ed_4yr,
            low_bmi, ovwt_bmi, obese_bmi,
            maternal_age, any_smoker, smokeSH, mean_cpss, mean_epsd,
            male, gest_age_w) %>%
  as.matrix()
head(W)
##      latina_re black_re other_re ed_no_hs ed_hs ed_aa ed_4yr low_bmi ovwt_bmi
## [1,]         1        0        0        0     0     1      0       0        0
## [2,]         0        0        1        0     0     1      0       0        0
## [3,]         0        0        0        0     0     0      0       0        0
## [4,]         0        0        0        0     0     1      0       0        0
## [5,]         0        1        0        0     0     0      1       0        0
## [6,]         1        0        0        0     0     1      0       0        0
##      obese_bmi maternal_age any_smoker smokeSH mean_cpss mean_epsd male
## [1,]         0           19          0       1        29         0    0
## [2,]         0           36          0       0        19         2    1
## [3,]         0           34          0       0        19         1    0
## [4,]         0           28          0       0        20         0    0
## [5,]         0           30          0       0        15         0    1
## [6,]         0           22          0       0        17         1    0
##      gest_age_w
## [1,]   40.57143
## [2,]   35.85714
## [3,]   40.42857
## [4,]   36.28571
## [5,]   38.42857
## [6,]   40.71429

Scaled the non-binary (continuous) covariates

W.s <- apply(W[,c(11, 14, 15, 17)], 2, scale) #' just the continuous ones

W.scaled <- cbind(W[,1:10], W.s[,1],
                  W[,12:13], W.s[,2:3],
                  W[,16], W.s[,4])
colnames(W.scaled)
##  [1] "latina_re"  "black_re"   "other_re"   "ed_no_hs"   "ed_hs"     
##  [6] "ed_aa"      "ed_4yr"     "low_bmi"    "ovwt_bmi"   "obese_bmi" 
## [11] ""           "any_smoker" "smokeSH"    "mean_cpss"  "mean_epsd" 
## [16] ""           ""
colnames(W.scaled) <- colnames(W)
head(W.scaled)
##      latina_re black_re other_re ed_no_hs ed_hs ed_aa ed_4yr low_bmi ovwt_bmi
## [1,]         1        0        0        0     0     1      0       0        0
## [2,]         0        0        1        0     0     1      0       0        0
## [3,]         0        0        0        0     0     0      0       0        0
## [4,]         0        0        0        0     0     1      0       0        0
## [5,]         0        1        0        0     0     0      1       0        0
## [6,]         1        0        0        0     0     1      0       0        0
##      obese_bmi maternal_age any_smoker smokeSH  mean_cpss  mean_epsd male
## [1,]         0  -1.39815187          0       1  3.3147856 -1.2832098    0
## [2,]         0   1.35109608          0       0  0.1179652 -0.6860171    1
## [3,]         0   1.02765515          0       0  0.1179652 -0.9846134    0
## [4,]         0   0.05733234          0       0  0.4376472 -1.2832098    0
## [5,]         0   0.38077328          0       0 -1.1607630 -1.2832098    1
## [6,]         0  -0.91299047          0       0 -0.5213989 -0.9846134    0
##      gest_age_w
## [1,]  0.7037686
## [2,] -1.9146645
## [3,]  0.6244221
## [4,] -1.6766251
## [5,] -0.4864283
## [6,]  0.7831150
summary(W.scaled)
##    latina_re         black_re         other_re          ed_no_hs     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :0.00000   Median :0.0000  
##  Mean   :0.2653   Mean   :0.1717   Mean   :0.06689   Mean   :0.1527  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000  
##      ed_hs            ed_aa            ed_4yr          low_bmi       
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000  
##  Median :0.0000   Median :0.0000   Median :0.0000   Median :0.00000  
##  Mean   :0.1851   Mean   :0.2319   Mean   :0.2185   Mean   :0.03344  
##  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
##     ovwt_bmi       obese_bmi       maternal_age        any_smoker     
##  Min.   :0.000   Min.   :0.0000   Min.   :-1.88331   Min.   :0.00000  
##  1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:-0.91299   1st Qu.:0.00000  
##  Median :0.000   Median :0.0000   Median : 0.05733   Median :0.00000  
##  Mean   :0.262   Mean   :0.1996   Mean   : 0.00000   Mean   :0.08696  
##  3rd Qu.:1.000   3rd Qu.:0.0000   3rd Qu.: 0.70421   3rd Qu.:0.00000  
##  Max.   :1.000   Max.   :1.0000   Max.   : 2.64486   Max.   :1.00000  
##     smokeSH         mean_cpss         mean_epsd            male       
##  Min.   :0.0000   Min.   :-5.9560   Min.   :-1.2832   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:-0.5214   1st Qu.:-0.7855   1st Qu.:0.0000  
##  Median :0.0000   Median : 0.0114   Median :-0.1884   Median :1.0000  
##  Mean   :0.2575   Mean   : 0.0000   Mean   : 0.0000   Mean   :0.5117  
##  3rd Qu.:1.0000   3rd Qu.: 0.5442   3rd Qu.: 0.6079   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   : 4.5935   Max.   : 6.0324   Max.   :1.0000  
##    gest_age_w     
##  Min.   :-7.7070  
##  1st Qu.:-0.3277  
##  Median : 0.1483  
##  Mean   : 0.0000  
##  3rd Qu.: 0.6244  
##  Max.   : 2.9255

Variance and histograms for the scaled covariates

var(W.scaled)
##                   latina_re      black_re      other_re     ed_no_hs
## latina_re     0.19514701784 -0.0456034002 -0.0177675585  0.039787884
## black_re     -0.04560340022  0.1423669175 -0.0114966555  0.012811803
## other_re     -0.01776755853 -0.0114966555  0.0624850693 -0.003531116
## ed_no_hs      0.03978788422  0.0128118032 -0.0035311156  0.129548893
## ed_hs         0.02896808807  0.0150675864  0.0065807155 -0.028296206
## ed_aa         0.01318258282  0.0192967133  0.0023291925 -0.035455487
## ed_4yr       -0.03571926262 -0.0163503842 -0.0034713927 -0.033409978
## low_bmi       0.00004479216  0.0020641722 -0.0011235368 -0.003997701
## ovwt_bmi      0.02081218148 -0.0003857103 -0.0052668120  0.004584976
## obese_bmi     0.02065416468  0.0069962872  0.0022620043  0.006318184
## maternal_age -0.10858833047 -0.0895968722 -0.0133074822 -0.137612465
## any_smoker   -0.00858889752  0.0185364907 -0.0013586957  0.017954193
## smokeSH      -0.00590510033  0.0349789477  0.0084246596  0.030936455
## mean_cpss    -0.04668314421 -0.0233123370  0.0255604593 -0.056115354
## mean_epsd     0.04322362524  0.0220578552  0.0202679881  0.067010728
## male         -0.00087717989  0.0002202281  0.0003322086 -0.013508570
## gest_age_w   -0.02639363168 -0.0368471269 -0.0054750903 -0.007692848
##                       ed_hs        ed_aa       ed_4yr        low_bmi
## latina_re     0.02896808807  0.013182583 -0.035719263  0.00004479216
## black_re      0.01506758640  0.019296713 -0.016350384  0.00206417224
## other_re      0.00658071548  0.002329193 -0.003471393 -0.00112353679
## ed_no_hs     -0.02829620561 -0.035455487 -0.033409978 -0.00399770067
## ed_hs         0.15098194378 -0.042960663 -0.040482163 -0.00173196369
## ed_aa        -0.04296066253  0.178312629 -0.050724638  0.00786102484
## ed_4yr       -0.04048216276 -0.050724638  0.170951784 -0.00173569637
## low_bmi      -0.00173196369  0.007861025 -0.001735696  0.03236233875
## ovwt_bmi     -0.00612657270  0.015075052  0.004074843 -0.00877179885
## obese_bmi     0.02441297380  0.009478520 -0.012402453 -0.00668149785
## maternal_age -0.10868310364 -0.040296710  0.109104452 -0.01089529356
## any_smoker    0.00174689441  0.011063665 -0.015673525  0.00155279503
## smokeSH       0.02483352246  0.022806677 -0.036244326  0.00142215122
## mean_cpss    -0.03433895561  0.030714793  0.028242932  0.00484169668
## mean_epsd     0.02014619096  0.025424770 -0.046368432  0.00946748435
## male          0.00006345557  0.000630823  0.006367953 -0.00262407430
## gest_age_w   -0.01723767893 -0.035168407  0.030147649 -0.00601412877
##                   ovwt_bmi    obese_bmi maternal_age   any_smoker      smokeSH
## latina_re     0.0208121815  0.020654165 -0.108588330 -0.008588898 -0.005905100
## black_re     -0.0003857103  0.006996287 -0.089596872  0.018536491  0.034978948
## other_re     -0.0052668120  0.002262004 -0.013307482 -0.001358696  0.008424660
## ed_no_hs      0.0045849757  0.006318184 -0.137612465  0.017954193  0.030936455
## ed_hs        -0.0061265727  0.024412974 -0.108683104  0.001746894  0.024833522
## ed_aa         0.0150750518  0.009478520 -0.040296710  0.011063665  0.022806677
## ed_4yr        0.0040748427 -0.012402453  0.109104452 -0.015673525 -0.036244326
## low_bmi      -0.0087717989 -0.006681498 -0.010895294  0.001552795  0.001422151
## ovwt_bmi      0.1935643614 -0.052338400  0.008900228 -0.006065606 -0.010623208
## obese_bmi    -0.0523383998  0.159910515  0.002790074  0.002717391  0.011052467
## maternal_age  0.0089002276  0.002790074  1.000000000 -0.046629611 -0.155964054
## any_smoker   -0.0060656056  0.002717391 -0.046629611  0.079483696  0.049010093
## smokeSH      -0.0106232083  0.011052467 -0.155964054  0.049010093  0.191419314
## mean_cpss    -0.0082476896 -0.008722611  0.100637638  0.017642908  0.031721118
## mean_epsd    -0.0019127169  0.028133038 -0.160410684  0.042144665  0.108180210
## male          0.0008361204 -0.001780489  0.023413804  0.002329193  0.002004449
## gest_age_w   -0.0148466585 -0.021461430  0.091663607 -0.014981418 -0.050311537
##                 mean_cpss    mean_epsd           male   gest_age_w
## latina_re    -0.046683144  0.043223625 -0.00087717989 -0.026393632
## black_re     -0.023312337  0.022057855  0.00022022814 -0.036847127
## other_re      0.025560459  0.020267988  0.00033220855 -0.005475090
## ed_no_hs     -0.056115354  0.067010728 -0.01350857023 -0.007692848
## ed_hs        -0.034338956  0.020146191  0.00006345557 -0.017237679
## ed_aa         0.030714793  0.025424770  0.00063082298 -0.035168407
## ed_4yr        0.028242932 -0.046368432  0.00636795270  0.030147649
## low_bmi       0.004841697  0.009467484 -0.00262407430 -0.006014129
## ovwt_bmi     -0.008247690 -0.001912717  0.00083612040 -0.014846658
## obese_bmi    -0.008722611  0.028133038 -0.00178048853 -0.021461430
## maternal_age  0.100637638 -0.160410684  0.02341380415  0.091663607
## any_smoker    0.017642908  0.042144665  0.00232919255 -0.014981418
## smokeSH       0.031721118  0.108180210  0.00200444935 -0.050311537
## mean_cpss     1.000000000  0.455187203 -0.00331530432 -0.037142336
## mean_epsd     0.455187203  1.000000000  0.00154181454 -0.137187808
## male         -0.003315304  0.001541815  0.25014184185 -0.007427180
## gest_age_w   -0.037142336 -0.137187808 -0.00742717951  1.000000000
ggplot(pivot_longer(as.data.frame(W.scaled), latina_re:gest_age_w, 
                    names_to = "exp", values_to = "value")) + 
    geom_histogram(aes(x = value)) + 
    facet_wrap(~ exp, scales = "free")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.3 Response data: birth weight

Y <- select(hs_data2, birth_weight) %>%
  as.matrix()
head(Y)
##      birth_weight
## [1,]         2860
## [2,]         2755
## [3,]         3505
## [4,]         2695
## [5,]         3355
## [6,]         3810

Distribution of birth weight and scaled birth weight

hist(Y, breaks = 20)

hist(scale(Y), breaks = 20)

1.4 Scatterplots of exposures and outcome (birth weight)

Both birth weight (Y) and the exposures are scaled here

NOTE: Don’t use these plots as a way to estimate how many predictors might make the cut. This should be done a priori

df <- as.data.frame(cbind(scale(Y), X.scaled))
# par(mfrow=c(5,4))
sapply(2:length(df), function(x){
  lm.x <- lm(birth_weight ~ df[,x], data = df)
  plot(df[,c(x, 1)],
       xlab = paste0(colnames(df)[x], " beta: ",
                     round(summary(lm.x)$coef[2,1],4),
                     "; p = ",
                     round(summary(lm.x)$coef[2,4],4)))
  abline(lm.x)
})

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
## 
## [[13]]
## NULL
## 
## [[14]]
## NULL
## 
## [[15]]
## NULL
## 
## [[16]]
## NULL
## 
## [[17]]
## NULL
## 
## [[18]]
## NULL
## 
## [[19]]
## NULL

2 Is gestational age a potenital mediator?

I.e., is there a relationship between our exposures and gestational age?

The DAG might look something like this:

exposures —> gestational age —> birth weight

2.1 Scatter plots for exposures and gestational age

Both gestational age and the exposures are scaled here. Gestational age measured in weeks from estimated date of conception to delivery

Since there were some (small) relationships between exposures and gestational age (based on simple linear regression models– namely the ozone and SES indicators), I’m going to omit this covariate for now.

df2 <- as.data.frame(cbind(W.scaled[,"gest_age_w"], X.scaled))
colnames(df2)[1] <- "gest_age_w"
# par(mfrow=c(5,4))
sapply(2:length(df2), function(x){
  lm.x <- lm(gest_age_w ~ df2[,x], data = df2)
  plot(df2[,c(x, 1)],
       xlab = paste0(colnames(df2)[x], " beta: ",
                     round(summary(lm.x)$coef[2,1],4),
                     "; p = ",
                     round(summary(lm.x)$coef[2,4],4)))
  abline(lm.x)
})

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL
## 
## [[9]]
## NULL
## 
## [[10]]
## NULL
## 
## [[11]]
## NULL
## 
## [[12]]
## NULL
## 
## [[13]]
## NULL
## 
## [[14]]
## NULL
## 
## [[15]]
## NULL
## 
## [[16]]
## NULL
## 
## [[17]]
## NULL
## 
## [[18]]
## NULL
## 
## [[19]]
## NULL

Dropping gest_age_w from the covariates

W.scaled2 <- W.scaled[,-c(ncol(W.scaled))]
colnames(W.scaled2)
##  [1] "latina_re"    "black_re"     "other_re"     "ed_no_hs"     "ed_hs"       
##  [6] "ed_aa"        "ed_4yr"       "low_bmi"      "ovwt_bmi"     "obese_bmi"   
## [11] "maternal_age" "any_smoker"   "smokeSH"      "mean_cpss"    "mean_epsd"   
## [16] "male"

3 RIDGE regression

To see if there might be something going on, Lauren suggested a ridge regression with a small penalty.

set.seed(123)

library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loaded glmnet 3.0-2
lambda_seq <- 10^seq(2, -2, by = -.1)

#' Best lambda from CV
ridge_cv <- cv.glmnet(X.scaled, scale(Y), alpha = 0, lambda = lambda_seq)
plot(ridge_cv)

best_lambda <- ridge_cv$lambda.min
best_lambda
## [1] 1.258925
#' Fit the model using the best_lambda
bw_ridge <- glmnet(X.scaled, scale(Y), alpha = 0, lambda = best_lambda)
summary(bw_ridge)
##           Length Class     Mode   
## a0         1     -none-    numeric
## beta      19     dgCMatrix S4     
## df         1     -none-    numeric
## dim        2     -none-    numeric
## lambda     1     -none-    numeric
## dev.ratio  1     -none-    numeric
## nulldev    1     -none-    numeric
## npasses    1     -none-    numeric
## jerr       1     -none-    numeric
## offset     1     -none-    logical
## call       5     -none-    call   
## nobs       1     -none-    numeric

Ridge regression coefficients

coef(bw_ridge)
## 20 x 1 sparse Matrix of class "dgCMatrix"
##                                           s0
## (Intercept)         -0.000000000000000214871
## mean_pm              0.009001292934592152947
## mean_o3             -0.048621510737574098748
## pct_tree_cover       0.001342271308724626074
## pct_impervious      -0.012521694754762259516
## mean_aadt_intensity -0.005312589097233479454
## dist_m_tri          -0.001659989297196438852
## dist_m_npl           0.001245550042200224651
## dist_m_waste_site    0.008057207316642686287
## dist_m_major_emit   -0.002818043274053380360
## dist_m_cafo         -0.003389166474688925200
## dist_m_mine_well    -0.009950869706335926587
## cvd_rate_adj        -0.016027916742713192721
## res_rate_adj        -0.009686648757473777238
## violent_crime_rate  -0.005457924031598470650
## property_crime_rate -0.000364961614232602028
## pct_less_hs         -0.021599485837649336911
## pct_unemp           -0.042979714186389718356
## pct_limited_eng     -0.010315534078569870563
## pct_hh_pov          -0.010927666957598947128

Ridge regression predictions

ridge_pred <- predict(bw_ridge, newx = X.scaled)
plot(scale(Y), ridge_pred)

actual <- scale(Y)
preds <- ridge_pred
rsq <- 1 - (sum((preds - actual) ^ 2))/(sum((actual - mean(actual)) ^ 2))

The R2 value for this model is 0.03. Based on these results, it doesn’t look like there’s much here.

4 Nonparametric Bayesian Shrinkage (NPB): Birth weight

Still, we wanted to try to fit the NPB model with these data.

4.1 Finding the NPB priors

Start with Lauren’s from the example in the vignette

In an email from April 29, Lauren provided me with some additional guidance on finding the NPB priors:

  • Keep alpha.pi and beta.pi set to 1, and then let a.phi1 take values 1, 10, and 100 and see how the results change.
  • Keep a.phi1 set to 1 (or 10 or 100), and mess with alpha.pi and beta.pi.
    • Run the following code: alpha.pi=1 beta.pi=1 plot(density(rbeta(10000, alpha.pi, beta.pi))) and then change alpha.pi and beta.pi and see how it changes the prior distribution.
    • This is the distribution of the probability of a main effect regression coefficient being 0 (aka exclusion probability). We don’t want this to be too informative (you don’t want high mass around just a few values). Also alpha.pi and beta.pi don’t have to be the same value. You might try alpha.pi = 1 and beta.pi = 2 to get a slightly lower prior probability of exclusion. Try not to change all three (alpha.pi, beta.pi, and a.phi1) at once.
  • When playing with the priors, set “interact=FALSE” and just fit the modelwith the main effects. Most of the interactions were null anyway so it shouldn’t change the results too much and it will make the code run a lot faster. Then when you find a set of priors you like, you can add in “interact=TRUE” and “XWinteract=TRUE.”

Some additional feedback from Lauren during our 6/10 meeting:

  • The confidence intervals were really wide and heavily skewed. I’m going to try adjusting the sig2inv.mu1 parameter after the a.phi1 parameter to see if this helps
  • the rbeta distributions is interpreted as the exclusion probability, so I should try to aim to have most of the mass of that distribution in the middle, since we hypothesize that maybe 40-60% of the predictors will be important. The way to do this is to set alpha.pi and beta.pi to the same value
  • I should set the burn in number to be about half the iterations

4.1.1 Vignette Priors

set.seed(123)

priors.npb.1 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1)

fit.npb.1 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.1, interact = F)
npb.sum.1 <- summary(fit.npb.1)
npb.sum.1$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.42293623  3.1955064     0.000000            0 0.022
##  [2,]    -7.77734425 17.0869661   -61.452320            0 0.208
##  [3,]    -0.34921645  2.9046941     0.000000            0 0.020
##  [4,]    -0.72527602  4.5329148    -9.113312            0 0.032
##  [5,]    -0.10902198  1.6231693     0.000000            0 0.008
##  [6,]    -0.22439232  3.4307474     0.000000            0 0.022
##  [7,]     0.13813603  3.2405735     0.000000            0 0.014
##  [8,]     0.50835031  5.5683575     0.000000            0 0.018
##  [9,]    -0.01225732  0.1949014     0.000000            0 0.004
## [10,]    -0.08255604  1.1888837     0.000000            0 0.006
## [11,]    -1.04104994  5.8717128   -21.539810            0 0.050
## [12,]    -0.44293941  3.6077235     0.000000            0 0.020
## [13,]    -1.34313104  6.6180314   -26.725959            0 0.050
## [14,]    -0.21174177  1.9565687     0.000000            0 0.016
## [15,]    -0.67619819  4.4923381    -7.070728            0 0.034
## [16,]    -0.69671736  4.5842905    -3.417542            0 0.026
## [17,]    -3.11238035 10.7411492   -42.631782            0 0.096
## [18,]    -0.58256651  4.1194858    -3.681873            0 0.032
## [19,]    -0.41969750  3.4549837     0.000000            0 0.020
plot(fit.npb.1$beta[,1], type = "l")

plot(fit.npb.1$beta[,2], type = "l")

plot(fit.npb.1$beta[,13], type = "l")

4.1.2 Adjust sig2inv.mu1

The default value for this parameter is 1. Adjusting this parameter may help narrow these confidence intervals a bit

4.1.2.1 Try sig2inv.mu1 = 10

That seemed to make the confidence intervals wider

priors.npb.2 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     sig2inv.mu1 = 10)

fit.npb.2 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.2, interact = F)
npb.sum.2 <- summary(fit.npb.2)
npb.sum.2$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.57320485  3.542448   -10.457164     0.000000 0.060
##  [2,]    -9.76615860 18.286535   -58.159264     0.000000 0.294
##  [3,]    -0.48674498  2.686788    -8.315099     0.000000 0.056
##  [4,]    -1.17031436  6.050804   -16.685042     0.000000 0.084
##  [5,]    -0.18426760  3.095979    -4.889139     0.000000 0.042
##  [6,]    -0.17137354  3.496180    -7.318933     0.000000 0.058
##  [7,]     0.12313460  3.434666    -1.841864     0.000000 0.054
##  [8,]     0.68066167  6.178553     0.000000     7.296187 0.050
##  [9,]     0.07551352  3.152292     0.000000     0.000000 0.038
## [10,]    -0.24556275  1.952749    -6.237070     0.000000 0.048
## [11,]    -1.22843970  5.154892   -17.980677     0.000000 0.096
## [12,]    -0.62684355  3.672637    -9.382015     0.000000 0.058
## [13,]    -2.42125423  8.397953   -34.673034     0.000000 0.122
## [14,]    -0.43994939  3.003203    -8.243475     0.000000 0.052
## [15,]    -1.12213230  4.786993   -16.736005     0.000000 0.078
## [16,]    -1.20895702  5.295603   -18.793484     0.000000 0.086
## [17,]    -4.07327965 11.899635   -49.929697     0.000000 0.158
## [18,]    -0.47703166  2.789470    -8.900533     0.000000 0.064
## [19,]    -0.81382810  4.258625   -11.526663     0.000000 0.070
plot(fit.npb.2$beta[,1], type = "l")

plot(fit.npb.2$beta[,2], type = "l")

plot(fit.npb.2$beta[,13], type = "l")

4.1.2.2 Try sig2inv.mu1 = 100

That also seemed to make the confidence intervals wider

priors.npb.3 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     sig2inv.mu1 = 100)

fit.npb.3 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.3, interact = F)
npb.sum.3 <- summary(fit.npb.3)
npb.sum.3$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.095131290  3.140375     0.000000    0.0000000 0.020
##  [2,]  -11.681887793 20.022155   -59.403940    0.0000000 0.292
##  [3,]    0.016920655  3.934123     0.000000    0.0000000 0.028
##  [4,]   -0.914938821  5.470093   -18.182994    0.0000000 0.040
##  [5,]   -0.248071932  2.846108     0.000000    0.0000000 0.034
##  [6,]   -0.059851371  2.346589     0.000000    0.0000000 0.030
##  [7,]    0.517100918  4.544226     0.000000    0.2186056 0.032
##  [8,]    1.371237081  8.156704     0.000000   17.0215314 0.050
##  [9,]    0.252087787  2.563337     0.000000    0.0000000 0.026
## [10,]   -0.003136574  2.221280     0.000000    0.0000000 0.030
## [11,]   -0.871617648  5.610285   -18.258161    0.0000000 0.040
## [12,]   -0.605872070  4.903910   -12.928430    0.0000000 0.044
## [13,]   -1.138512749  6.681420   -23.502834    0.0000000 0.058
## [14,]   -0.053623560  2.120945     0.000000    0.0000000 0.024
## [15,]   -0.430496025  4.375130    -2.256698    0.0000000 0.042
## [16,]   -0.556105327  7.871809   -22.192950    0.0000000 0.058
## [17,]   -4.035417542 12.806278   -49.923011    0.0000000 0.120
## [18,]   -0.334067046  3.709821     0.000000    0.0000000 0.036
## [19,]   -1.149915311  7.331426   -10.142257    0.0000000 0.030
plot(fit.npb.3$beta[,1], type = "l")

plot(fit.npb.3$beta[,2], type = "l")

plot(fit.npb.3$beta[,13], type = "l")

4.1.2.3 Try sig2inv.mu1 = 0.1?

This might be too narrow?

priors.npb.4 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     sig2inv.mu1 = 0.1)

fit.npb.4 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.4, interact = F)
npb.sum.4 <- summary(fit.npb.4)
npb.sum.4$main.effects
##       Posterior Mean          SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.105976545  1.92589731      0.00000            0 0.008
##  [2,]   -3.708253022 12.88063652    -51.92299            0 0.086
##  [3,]   -0.004328270  0.09383458      0.00000            0 0.004
##  [4,]   -0.133231949  2.20638272      0.00000            0 0.006
##  [5,]   -0.010487186  0.16648611      0.00000            0 0.006
##  [6,]   -0.003353958  0.07499678      0.00000            0 0.002
##  [7,]    0.000000000  0.00000000      0.00000            0 0.000
##  [8,]    0.226086189  3.07843115      0.00000            0 0.010
##  [9,]   -0.015321283  1.31966148      0.00000            0 0.006
## [10,]   -0.127552312  1.86015024      0.00000            0 0.008
## [11,]   -0.385794753  4.12927001      0.00000            0 0.010
## [12,]   -0.400850294  4.05400654      0.00000            0 0.014
## [13,]   -0.786528947  5.60881484      0.00000            0 0.024
## [14,]   -0.112417701  1.86716884      0.00000            0 0.004
## [15,]   -0.266540085  3.04248660      0.00000            0 0.012
## [16,]   -0.063191855  1.00336679      0.00000            0 0.004
## [17,]   -1.612479624  8.44032992    -36.35425            0 0.040
## [18,]   -0.098734619  1.42912623      0.00000            0 0.008
## [19,]   -0.238959635  3.10000213      0.00000            0 0.008
plot(fit.npb.4$beta[,1], type = "l")

plot(fit.npb.4$beta[,2], type = "l")

plot(fit.npb.4$beta[,13], type = "l")

4.1.2.4 Try sig2inv.mu1 = 0.5?

Also too narrow– going back to 1

priors.npb.4 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     sig2inv.mu1 = 0.1)

fit.npb.4 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.4, interact = F)
npb.sum.4 <- summary(fit.npb.4)
npb.sum.4$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.16653792  2.1552452     0.000000            0 0.016
##  [2,]    -3.07753903 11.8843477   -50.000667            0 0.082
##  [3,]    -0.05255064  0.9993058     0.000000            0 0.006
##  [4,]    -0.28271383  2.5061098     0.000000            0 0.022
##  [5,]    -0.06295496  1.1247119     0.000000            0 0.006
##  [6,]     0.02999439  0.7687866     0.000000            0 0.006
##  [7,]     0.21649613  3.4418344     0.000000            0 0.008
##  [8,]     0.50724409  4.6522999     0.000000            0 0.014
##  [9,]    -0.03404257  0.5868220     0.000000            0 0.004
## [10,]    -0.12325749  1.8669486     0.000000            0 0.014
## [11,]    -0.20457023  2.3663177     0.000000            0 0.010
## [12,]    -0.29222742  3.2766185     0.000000            0 0.016
## [13,]    -0.41257543  3.4489993     0.000000            0 0.020
## [14,]    -0.13817717  1.3110905     0.000000            0 0.012
## [15,]    -0.21872498  3.3383078     0.000000            0 0.020
## [16,]    -0.24942108  2.5110565     0.000000            0 0.014
## [17,]    -0.64575618  4.5204832    -4.510545            0 0.026
## [18,]    -0.25473085  2.7429636     0.000000            0 0.016
## [19,]    -0.40524422  2.9289043     0.000000            0 0.024
plot(fit.npb.4$beta[,1], type = "l")

plot(fit.npb.4$beta[,2], type = "l")

plot(fit.npb.4$beta[,13], type = "l")

4.1.3 Next, try adjusting a.phi1

This parameter is the shape parameter for gamma prior on phi sub1 inv2. It’s default value is 1

4.1.3.1 Try making a.phi1 = 10

PIPs again increase slightly

priors.npb.5 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10)

fit.npb.5 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.5, interact = F)
npb.sum.5 <- summary(fit.npb.5)
npb.sum.5$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.581552391  3.821619   -6.9203902      0.00000 0.040
##  [2,]  -11.473528795 20.634360  -61.5619258      0.00000 0.282
##  [3,]   -0.280908392  2.581600    0.0000000      0.00000 0.026
##  [4,]   -0.685781042  4.763465  -11.8478088      0.00000 0.048
##  [5,]    0.007128347  1.470303    0.0000000      0.00000 0.020
##  [6,]   -0.063487381  2.012826    0.0000000      0.00000 0.026
##  [7,]    0.346319628  4.162866    0.0000000      0.00000 0.030
##  [8,]    1.383203139  7.315417    0.0000000     28.75221 0.050
##  [9,]    0.046533892  1.936927    0.0000000      0.00000 0.032
## [10,]   -0.223120797  2.664584    0.0000000      0.00000 0.024
## [11,]   -1.097947666  5.698955  -18.5451922      0.00000 0.056
## [12,]   -0.584701314  5.176760   -5.6821139      0.00000 0.036
## [13,]   -1.339277024  6.544944  -26.5626778      0.00000 0.058
## [14,]   -0.093226827  2.408309    0.0000000      0.00000 0.030
## [15,]   -0.525808339  3.716516   -5.3152818      0.00000 0.032
## [16,]   -0.851948191  5.505950  -13.1334253      0.00000 0.040
## [17,]   -4.359826211 13.001580  -50.7884608      0.00000 0.136
## [18,]   -0.251327971  2.127834   -0.5059922      0.00000 0.032
## [19,]   -0.965748807  6.050453  -14.9412688      0.00000 0.052
plot(fit.npb.5$beta[,1], type = "l")

plot(fit.npb.5$beta[,2], type = "l")

plot(fit.npb.5$beta[,13], type = "l")

4.1.3.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10

Definitely worse in terms of the confidence intervals

priors.npb.6 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 10)

fit.npb.6 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.6, interact = F)
npb.sum.6 <- summary(fit.npb.6)
npb.sum.6$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -0.4928015  3.286066    -4.847052      0.00000 0.034
##  [2,]    -17.4623243 22.302478   -64.918518      0.00000 0.452
##  [3,]     -0.2814630  2.585887    -2.677709      0.00000 0.034
##  [4,]     -2.0754808  8.100446   -33.727567      0.00000 0.088
##  [5,]      0.1666961  3.909889     0.000000      0.00000 0.032
##  [6,]     -0.0959787  3.649590     0.000000      0.00000 0.024
##  [7,]      0.2210554  3.615359     0.000000      0.00000 0.040
##  [8,]      2.2210034 10.266335     0.000000     40.50969 0.072
##  [9,]      0.1968973  2.996566     0.000000      0.00000 0.030
## [10,]     -0.1913085  2.891537    -4.718272      0.00000 0.038
## [11,]     -1.8156249  7.942042   -29.464176      0.00000 0.094
## [12,]     -0.7893540  4.525473   -12.672066      0.00000 0.052
## [13,]     -2.9937749 10.303755   -43.520723      0.00000 0.110
## [14,]     -0.2091624  3.483150    -8.291385      0.00000 0.050
## [15,]     -1.1753419  5.587383   -17.375845      0.00000 0.068
## [16,]     -1.5192580  7.132802   -26.391704      0.00000 0.070
## [17,]     -7.1958362 15.741948   -53.830439      0.00000 0.218
## [18,]     -0.3065605  4.195224    -8.740062      0.00000 0.046
## [19,]     -1.0956311  6.006086   -19.708031      0.00000 0.052
plot(fit.npb.6$beta[,1], type = "l")

plot(fit.npb.6$beta[,2], type = "l")

plot(fit.npb.6$beta[,13], type = "l")

4.1.3.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100

Not an improvement

priors.npb.7 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 100)

fit.npb.7 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.7, interact = F)
npb.sum.7 <- summary(fit.npb.7)
npb.sum.7$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.22104871  2.757115     0.000000      0.00000 0.038
##  [2,]   -20.41085141 23.151185   -66.092224      0.00000 0.524
##  [3,]    -0.24631131  2.661488     0.000000      0.00000 0.034
##  [4,]    -2.09887813  8.881620   -35.022901      0.00000 0.072
##  [5,]     0.11682112  4.046631     0.000000      0.00000 0.032
##  [6,]     0.20596362  2.984464     0.000000      0.00000 0.028
##  [7,]     0.38192374  3.429390     0.000000      0.00000 0.038
##  [8,]     2.55506778 10.355460     0.000000     44.07723 0.082
##  [9,]     0.26757819  3.207645     0.000000      0.00000 0.030
## [10,]    -0.21437574  3.450039     0.000000      0.00000 0.036
## [11,]    -1.88737645  7.711836   -30.378845      0.00000 0.074
## [12,]    -1.04124580  7.530258   -21.971510      0.00000 0.084
## [13,]    -3.43025706 11.053184   -43.266666      0.00000 0.122
## [14,]     0.03994552  4.879433    -5.159376      0.00000 0.046
## [15,]    -1.57907973  7.007126   -23.564577      0.00000 0.084
## [16,]    -1.07553220  6.226368   -21.246083      0.00000 0.064
## [17,]    -7.66359663 16.478164   -57.107719      0.00000 0.218
## [18,]    -0.62459729  4.766040   -14.155347      0.00000 0.050
## [19,]    -0.94616666  5.567583   -17.775554      0.00000 0.058
plot(fit.npb.7$beta[,1], type = "l")

plot(fit.npb.7$beta[,2], type = "l")

plot(fit.npb.7$beta[,13], type = "l")

4.1.3.4 Try making a.phi1 = 10 and sig2inv.mu1 = 0.1

Nope

priors.npb.8 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 0.1)

fit.npb.8 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.8, interact = F)
npb.sum.8 <- summary(fit.npb.8)
npb.sum.8$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.15754855  2.1361739     0.000000            0 0.016
##  [2,]    -6.30991048 16.1308085   -56.736167            0 0.160
##  [3,]    -0.11841473  1.8581095     0.000000            0 0.018
##  [4,]    -0.34787660  3.0322241     0.000000            0 0.020
##  [5,]    -0.04711339  0.7135716     0.000000            0 0.008
##  [6,]     0.24920300  3.4996638     0.000000            0 0.024
##  [7,]     0.01351782  0.7029928     0.000000            0 0.010
##  [8,]     0.39877502  3.9120930     0.000000            0 0.028
##  [9,]     0.11288992  2.0333384     0.000000            0 0.020
## [10,]    -0.09659308  1.3557107     0.000000            0 0.014
## [11,]    -0.55699468  4.2665677     0.000000            0 0.026
## [12,]    -0.74334034  5.3325922    -9.749636            0 0.044
## [13,]    -1.22844989  6.3929916   -23.584655            0 0.044
## [14,]    -0.01716576  1.0229667     0.000000            0 0.014
## [15,]    -0.35293666  3.0341882     0.000000            0 0.018
## [16,]    -0.50600924  3.9284400     0.000000            0 0.030
## [17,]    -1.72631176  8.0307307   -34.544925            0 0.058
## [18,]    -0.20194757  2.5711228     0.000000            0 0.020
## [19,]    -0.27606966  2.9828552     0.000000            0 0.016
plot(fit.npb.8$beta[,1], type = "l")

plot(fit.npb.8$beta[,2], type = "l")

plot(fit.npb.8$beta[,13], type = "l")

4.1.3.5 Try making a.phi1 = 100

priors.npb.9 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 1)

fit.npb.9 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.9, interact = F)
npb.sum.9 <- summary(fit.npb.9)
npb.sum.9$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.49081449  4.629401   -11.061650    0.0000000 0.072
##  [2,]   -14.61898497 21.824480   -64.166897    0.0000000 0.384
##  [3,]    -0.42735467  3.520918    -9.044472    0.0000000 0.080
##  [4,]    -1.32439326  7.050419   -24.968932    0.0000000 0.090
##  [5,]    -0.16896624  3.790224    -5.342744    0.0000000 0.064
##  [6,]     0.08680043  3.475061    -1.363339    0.7361937 0.054
##  [7,]     0.21923467  4.910740    -5.391763    7.9913316 0.074
##  [8,]     3.71035015 12.113573     0.000000   47.9510551 0.142
##  [9,]     0.80244892  5.594285     0.000000   15.8078721 0.064
## [10,]    -0.01292279  3.611167    -4.312884    0.9968086 0.064
## [11,]    -1.62260660  7.209524   -26.960978    0.0000000 0.094
## [12,]    -0.50405112  4.702056   -10.830387    0.2037986 0.080
## [13,]    -2.37992543  9.107271   -34.744576    0.0000000 0.108
## [14,]     0.14853855  4.109383    -2.460935    0.2037986 0.062
## [15,]    -0.61961911  3.762306    -8.832907    0.0000000 0.062
## [16,]    -0.96949695  6.278286   -20.260757    0.0000000 0.092
## [17,]    -5.73190550 14.774390   -52.281428    0.0000000 0.188
## [18,]    -0.82480285  4.947323   -17.562099    0.0000000 0.092
## [19,]    -1.39328755  7.002128   -22.543324    0.0000000 0.092
plot(fit.npb.9$beta[,1], type = "l")

plot(fit.npb.9$beta[,2], type = "l")

plot(fit.npb.9$beta[,13], type = "l")

4.1.3.6 Try making a.phi1 = 100 and sig2inv.mu1 = 10

priors.npb.10 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 10)

fit.npb.10 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.10, interact = F)
npb.sum.10 <- summary(fit.npb.10)
npb.sum.10$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -0.8955066  4.608069   -14.166614    2.0820725 0.156
##  [2,]    -22.3467639 22.833264   -67.216503    0.0000000 0.632
##  [3,]     -1.0899915  6.266198   -17.583821    4.5550541 0.162
##  [4,]     -2.3737495  8.805426   -33.820739    4.7554120 0.202
##  [5,]     -0.3559511  4.337999   -11.343352    6.0657719 0.142
##  [6,]      0.1227590  5.098899   -11.856142   13.5502226 0.150
##  [7,]      1.0923054  6.925488    -7.915575   22.8299239 0.154
##  [8,]      4.2048826 13.169188    -6.612679   45.9947819 0.230
##  [9,]      0.3450384  4.113284    -7.279409   11.8604506 0.138
## [10,]     -0.3254172  4.709016   -10.877562    7.9023656 0.140
## [11,]     -3.5999901  9.997178   -34.815423    0.6276695 0.224
## [12,]     -1.5518494  8.235853   -23.908544    6.9765096 0.192
## [13,]     -3.6389730 10.311734   -37.952368    2.9627703 0.252
## [14,]     -0.2175783  5.652138   -13.317473   11.8566847 0.158
## [15,]     -1.6873238  6.697213   -22.019093    1.7058149 0.162
## [16,]     -2.5428487  9.312482   -32.166143    6.5707112 0.210
## [17,]     -8.0393071 16.401686   -53.826942    0.0000000 0.310
## [18,]     -1.1660232  6.244023   -20.317334    6.9765096 0.176
## [19,]     -1.7402131  7.377813   -26.084992    3.0057293 0.182
plot(fit.npb.10$beta[,1], type = "l")

plot(fit.npb.10$beta[,2], type = "l")

plot(fit.npb.10$beta[,13], type = "l")

4.1.3.7 Try making a.phi1 = 100 and sig2inv.mu1 = 100

Still not great.

priors.npb.11 <- list(alpha.pi = 1, beta.pi = 1, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 100)

fit.npb.11 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.11, interact = F)
npb.sum.11 <- summary(fit.npb.11)
npb.sum.11$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -0.7950218  5.244945   -16.970623     4.537841 0.126
##  [2,]    -27.2058531 21.929493   -68.156952     0.000000 0.728
##  [3,]     -0.4288572  4.965463   -13.452307     3.912050 0.116
##  [4,]     -3.2572407  9.847259   -36.272125     1.027606 0.196
##  [5,]      0.2192220  4.831992    -5.568812    12.263959 0.132
##  [6,]      0.2665736  6.216597    -9.888779    14.244300 0.128
##  [7,]      1.0229973  6.669657    -4.270251    21.779034 0.148
##  [8,]      3.9530248 11.917094    -2.002081    44.514415 0.198
##  [9,]      1.2600757  7.541149    -8.009803    27.325068 0.172
## [10,]     -0.1415800  6.526332   -12.783186     8.690473 0.134
## [11,]     -3.3426308 10.322144   -38.047100     0.000000 0.172
## [12,]     -1.2593536  8.369107   -24.126005     9.884446 0.164
## [13,]     -3.5622279 11.517625   -40.880095     5.691475 0.212
## [14,]      0.5091027  6.421518   -13.023474    23.268069 0.156
## [15,]     -2.4704408  8.688523   -31.742324     4.421180 0.198
## [16,]     -2.6022436 10.254831   -36.585588     5.991238 0.182
## [17,]    -10.5991456 17.509123   -51.783170     0.000000 0.362
## [18,]     -1.4959087  7.371900   -28.144047     4.255691 0.142
## [19,]     -2.1099624  8.284785   -28.971796     3.879175 0.160
plot(fit.npb.11$beta[,1], type = "l")

plot(fit.npb.11$beta[,2], type = "l")

plot(fit.npb.11$beta[,13], type = "l")

4.1.4 Adjust alpha.pi and beta.pi

For now, leave a.phi1 and sig2inv.mu1 alone for now.

alpha.pi and beta.pi are responisble for the exclusion probability distribution. If we thing we want ~50% of our covariates, we need the mass of this distribution to be somewhere between 0.4 and 0.6. To do this, we set alpha.pi and beta.pi to the same value

4.1.4.1 Try making alpha.pi and beta.pi 2

plot(density(rbeta(10000, 2, 2)))

priors.npb.12 <- list(alpha.pi = 2, beta.pi = 2, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 1)

fit.npb.12 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.12, interact = F)
npb.sum.12 <- summary(fit.npb.12)
npb.sum.12$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.548121233  4.0985444    -8.338155      0.00000 0.040
##  [2,]  -11.687879813 19.4861619   -59.865735      0.00000 0.330
##  [3,]   -0.377632766  3.5324695    -1.735183      0.00000 0.034
##  [4,]   -1.194321283  5.7745675   -22.468250      0.00000 0.060
##  [5,]   -0.054825515  0.8274353     0.000000      0.00000 0.010
##  [6,]    0.062375455  2.4983039     0.000000      0.00000 0.028
##  [7,]   -0.036286670  2.1712654     0.000000      0.00000 0.024
##  [8,]    1.280627308  7.7689013     0.000000     31.03455 0.048
##  [9,]    0.008804336  0.6536777     0.000000      0.00000 0.008
## [10,]   -0.196413518  1.6544415     0.000000      0.00000 0.020
## [11,]   -1.760859372  7.2421988   -28.993435      0.00000 0.072
## [12,]   -0.836066000  4.7927503   -16.370721      0.00000 0.056
## [13,]   -1.918927858  7.8922310   -30.315590      0.00000 0.084
## [14,]   -0.316961255  2.2995747     0.000000      0.00000 0.024
## [15,]   -0.842456591  5.1178116   -13.396083      0.00000 0.038
## [16,]   -1.653424484  6.8653248   -27.957693      0.00000 0.070
## [17,]   -5.446856463 13.2720490   -46.601885      0.00000 0.180
## [18,]   -0.512403416  3.4378903    -6.786906      0.00000 0.038
## [19,]   -0.772038898  4.4605108   -13.288182      0.00000 0.044
plot(fit.npb.12$beta[,1], type = "l")

plot(fit.npb.12$beta[,2], type = "l")

plot(fit.npb.12$beta[,13], type = "l")

4.1.4.2 Try making alpha.pi and beta.pi 5

plot(density(rbeta(10000, 5, 5)))

priors.npb.13 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 1)

fit.npb.13 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.13, interact = F)
npb.sum.13 <- summary(fit.npb.13)
npb.sum.13$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.5175818  6.596365    -24.41351    0.0000000 0.224
##  [2,]    -17.5507282 19.138119    -62.21890    0.0000000 0.658
##  [3,]     -1.3631356  4.765574    -16.61580    0.0000000 0.172
##  [4,]     -2.8709572  7.662579    -26.46499    0.0000000 0.234
##  [5,]     -1.0397884  4.671120    -15.15083    0.6703999 0.138
##  [6,]     -0.9327580  5.341849    -16.01600    3.3862704 0.168
##  [7,]     -0.4536953  4.155734    -11.29585    1.9789806 0.112
##  [8,]      0.7424446  8.219155    -10.60626   28.5344189 0.136
##  [9,]     -0.3941859  4.179095    -11.10859    4.1061592 0.132
## [10,]     -0.7716985  4.078803    -12.39498    2.0461425 0.128
## [11,]     -4.6362502 10.539032    -39.73617    0.4258184 0.284
## [12,]     -2.4914016  7.135969    -27.08012    0.0000000 0.224
## [13,]     -4.9526503 10.546835    -38.60816    0.0000000 0.302
## [14,]     -0.9026749  4.814544    -14.34211    1.3585607 0.150
## [15,]     -3.1748864  8.047856    -31.50061    0.0000000 0.238
## [16,]     -2.8539683  7.600540    -25.40936    1.1146465 0.252
## [17,]     -9.3501374 15.413350    -51.29416    0.0000000 0.420
## [18,]     -1.8299513  5.819704    -18.76603    0.0000000 0.188
## [19,]     -2.3426401  6.320223    -24.31243    0.0000000 0.218
plot(fit.npb.13$beta[,1], type = "l")

plot(fit.npb.13$beta[,2], type = "l")

plot(fit.npb.13$beta[,13], type = "l")

4.1.4.3 Try making alpha.pi and beta.pi 8

plot(density(rbeta(10000, 8, 8)))

priors.npb.14 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 1)

fit.npb.14 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.14, interact = F)
npb.sum.14 <- summary(fit.npb.14)
npb.sum.14$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.2208720  6.711495   -22.803458     1.487951 0.234
##  [2,]    -21.1917922 20.300690   -64.503294     0.000000 0.748
##  [3,]     -1.5420847  6.081411   -22.659259     2.424170 0.200
##  [4,]     -3.3387828  9.164850   -32.417900     0.000000 0.254
##  [5,]     -1.2164232  6.085029   -18.771920     6.704256 0.202
##  [6,]     -0.3928061  4.866547   -12.933445     8.727297 0.172
##  [7,]     -0.1886759  4.562750   -10.614897     9.006422 0.186
##  [8,]      1.9671743  9.803766    -9.345317    38.040782 0.210
##  [9,]      0.3059293  5.512448   -10.409413    13.616308 0.188
## [10,]     -1.0485428  5.747603   -16.448024     6.972009 0.210
## [11,]     -4.8200773 10.904677   -38.663265     3.040430 0.326
## [12,]     -2.2364080  8.259488   -24.910237     5.348292 0.268
## [13,]     -5.2746140 12.231653   -42.316979     1.903118 0.320
## [14,]     -0.9090102  4.871460   -13.365501     2.922642 0.174
## [15,]     -3.5391613  8.769814   -31.460777     2.242225 0.292
## [16,]     -3.4732758  9.609786   -33.340210     5.511336 0.310
## [17,]    -10.7033356 17.707997   -60.140820     0.000000 0.446
## [18,]     -2.3657150  7.688985   -26.275109     3.211347 0.240
## [19,]     -2.4531044  8.486316   -25.800574     5.974371 0.238
plot(fit.npb.14$beta[,1], type = "l")

plot(fit.npb.14$beta[,2], type = "l")

plot(fit.npb.14$beta[,13], type = "l")

4.1.4.4 Try making alpha.pi and beta.pi 10

This might be getting too big… Is there enough space in the tails?

plot(density(rbeta(10000, 10, 10)))

priors.npb.15 <- list(alpha.pi = 10, beta.pi = 10, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 1)

fit.npb.15 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.15, interact = F)
npb.sum.15 <- summary(fit.npb.15)
npb.sum.15$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.8857311  7.232977    -23.56255    5.3865842 0.306
##  [2,]    -17.1289824 18.507669    -63.63778    0.0000000 0.712
##  [3,]     -1.6323034  6.176281    -16.91383    6.2286586 0.238
##  [4,]     -3.6020419  8.267143    -29.03935    2.8766411 0.332
##  [5,]     -1.5654875  5.200797    -16.56090    3.6318501 0.214
##  [6,]     -0.6052462  7.281443    -17.65679   16.2820039 0.224
##  [7,]     -0.2546045  5.397410    -13.79061   12.1249885 0.198
##  [8,]      1.5759810  9.269974    -12.00373   37.5717956 0.224
##  [9,]     -0.3945213  6.360807    -14.74431   11.4964058 0.214
## [10,]     -1.2363866  5.434546    -16.68779    5.3865842 0.212
## [11,]     -5.1448347  9.714039    -31.62555    0.8892998 0.362
## [12,]     -3.1353175  7.871987    -26.39008    4.2840909 0.296
## [13,]     -4.7910311  9.203093    -29.48755    0.7573745 0.382
## [14,]     -0.9290989  6.399126    -14.91728    9.2503359 0.246
## [15,]     -3.1021597  7.106243    -23.80494    0.6839960 0.288
## [16,]     -4.2095545  8.910515    -30.64462    3.1182682 0.354
## [17,]     -9.8871161 14.791150    -51.45763    0.0000000 0.516
## [18,]     -2.5865340  7.240974    -25.16143    2.8766411 0.288
## [19,]     -2.6717984  8.003829    -24.29969    3.6711970 0.280
plot(fit.npb.15$beta[,1], type = "l")

plot(fit.npb.15$beta[,2], type = "l")

plot(fit.npb.15$beta[,13], type = "l")

4.1.5 Set alpha.pi and beta.pi to 8, readjust a.phi1 and sig2inv.mu1

Set alpha.pi and beta.pi to 8, try adjusting a.phi1 and sig2inv.mu1

4.1.5.1 Try making a.phi1 = 10 and sig2inv.mu1 = 1

priors.npb.16 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 1)

fit.npb.16 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.16, interact = F)
npb.sum.16 <- summary(fit.npb.16)
npb.sum.16$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.5584270  7.366529   -23.113622     3.844689 0.274
##  [2,]    -25.6541530 20.781097   -70.525111     0.000000 0.822
##  [3,]     -0.8456003  5.351849   -16.155950     6.628430 0.186
##  [4,]     -3.9716715  9.687439   -33.874152     3.100136 0.316
##  [5,]     -0.6566620  5.421212   -15.942915     6.547565 0.202
##  [6,]     -0.5472398  5.943481   -13.675801    12.959635 0.222
##  [7,]      0.5289133  7.935148   -12.285135    19.441791 0.220
##  [8,]      3.6669225 11.767144    -6.936092    44.027025 0.250
##  [9,]      0.9840048  7.205198   -11.178864    22.955402 0.224
## [10,]     -0.8924849  6.122219   -16.471644     9.911731 0.208
## [11,]     -5.2796511 11.725045   -39.869035     2.702883 0.346
## [12,]     -2.9070741  9.647720   -31.198810     8.953711 0.298
## [13,]     -5.6108588 12.693757   -42.736833     6.458987 0.360
## [14,]      0.1968864  9.101304   -16.269186    21.797660 0.232
## [15,]     -3.9612949 10.534515   -38.442630     3.788611 0.298
## [16,]     -2.7523511  8.514096   -28.261148     5.919833 0.280
## [17,]    -13.8663411 18.848322   -57.570802     0.123517 0.534
## [18,]     -1.3469540  6.406979   -17.682461     9.557808 0.250
## [19,]     -2.8201205  9.091258   -30.279149     5.574324 0.270
plot(fit.npb.16$beta[,1], type = "l")

plot(fit.npb.16$beta[,2], type = "l")

plot(fit.npb.16$beta[,13], type = "l")

Same priors, but now with scaleY = F

priors.npb.16a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 1)

fit.npb.16a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.16a, interact = F)
npb.sum.16a <- summary(fit.npb.16a)
npb.sum.16a$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,] -0.01991776141 0.2863296   -0.6992964    0.5727729 0.464
##  [2,]  0.00056061079 0.3187154   -0.7958305    0.7398862 0.506
##  [3,] -0.03471161517 0.3055429   -0.7110587    0.6485434 0.480
##  [4,] -0.02622573229 0.2940953   -0.6886811    0.6722611 0.494
##  [5,] -0.03089337315 0.3053399   -0.7938389    0.6140577 0.486
##  [6,] -0.01576318329 0.3080216   -0.7080183    0.6828620 0.494
##  [7,] -0.03381090913 0.3125751   -0.8258076    0.6287432 0.478
##  [8,] -0.00861707349 0.3123751   -0.7661726    0.6777252 0.518
##  [9,] -0.02157509124 0.3382525   -0.8302362    0.6707648 0.546
## [10,] -0.02870830928 0.3003430   -0.7902254    0.6221704 0.524
## [11,] -0.02532884091 0.3069601   -0.7470368    0.6530465 0.502
## [12,] -0.02095698714 0.3255910   -0.7847839    0.6771903 0.498
## [13,] -0.02059436427 0.2980543   -0.7584723    0.6586088 0.502
## [14,] -0.01576474199 0.3142070   -0.6979910    0.6499558 0.510
## [15,] -0.02125723989 0.3174181   -0.7684473    0.6804516 0.488
## [16,]  0.00006688786 0.3072310   -0.7088393    0.7004202 0.486
## [17,] -0.02806268605 0.3062034   -0.7106627    0.6624262 0.488
## [18,] -0.00615786566 0.3145777   -0.7314019    0.6625514 0.492
## [19,] -0.02278172580 0.3291935   -0.7219907    0.7325880 0.512
plot(fit.npb.16a$beta[,1], type = "l")

plot(fit.npb.16a$beta[,2], type = "l")

plot(fit.npb.16a$beta[,13], type = "l")

4.1.5.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10

priors.npb.17 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 10)

fit.npb.17 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.17, interact = F)
npb.sum.17 <- summary(fit.npb.17)
npb.sum.17$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.1114031  7.367810   -25.931299     9.363107 0.274
##  [2,]    -25.3346069 19.863539   -66.334946     0.000000 0.848
##  [3,]     -1.2404926  5.768699   -17.263289     9.161327 0.256
##  [4,]     -3.8940779  9.888461   -32.943501     5.677994 0.312
##  [5,]     -0.6420384  5.274159   -13.895186     9.204381 0.230
##  [6,]     -0.6245669  6.644525   -17.225677    13.345719 0.256
##  [7,]      0.3270219  7.875332   -13.175990    17.337589 0.250
##  [8,]      4.2127861 13.452337   -10.933251    47.477502 0.306
##  [9,]      1.1309278  7.689213    -9.963829    23.003936 0.242
## [10,]     -1.5547519  7.118842   -19.923356    11.966823 0.242
## [11,]     -6.8228893 13.060407   -42.538875     5.610858 0.432
## [12,]     -2.2978183  9.767088   -28.140006    11.749363 0.304
## [13,]     -5.9216198 12.252662   -41.722235     6.660095 0.394
## [14,]     -0.3628362  6.824363   -15.578790    15.121946 0.240
## [15,]     -3.8462739  9.954491   -31.914553     6.102279 0.336
## [16,]     -3.0798149  9.305046   -29.956764    10.470342 0.340
## [17,]    -14.8435061 17.959995   -56.317799     0.000000 0.592
## [18,]     -1.5879665  7.736488   -21.215828    11.581411 0.268
## [19,]     -2.7080880  8.515935   -27.715093     8.802932 0.302
plot(fit.npb.17$beta[,1], type = "l")

plot(fit.npb.17$beta[,2], type = "l")

plot(fit.npb.17$beta[,13], type = "l")

4.1.5.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100

priors.npb.18 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 100)

fit.npb.18 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.18, interact = F)
npb.sum.18 <- summary(fit.npb.18)
npb.sum.18$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.7940023  7.665154   -25.139592     1.967637 0.262
##  [2,]    -24.6794003 19.177369   -64.105042     0.000000 0.826
##  [3,]     -2.1369699  7.197838   -23.287142     2.457191 0.234
##  [4,]     -4.0511400  9.299041   -31.811440     2.462465 0.330
##  [5,]     -0.6707083  5.768537   -16.919951     7.073050 0.168
##  [6,]     -0.6148905  6.696902   -17.526985    15.081068 0.212
##  [7,]     -0.1515304  8.041328   -15.751521    17.158266 0.226
##  [8,]      5.4699606 14.981968    -5.774078    51.032414 0.246
##  [9,]      0.6457981  7.848650   -11.650429    30.095922 0.180
## [10,]     -1.4251170  7.085737   -22.208937     7.116970 0.228
## [11,]     -6.3271767 12.483270   -43.419739     2.008763 0.380
## [12,]     -3.3099004 10.237558   -32.374922     4.621957 0.266
## [13,]     -5.9102346 11.934866   -39.143598     1.522695 0.350
## [14,]      0.1825934 11.150450   -19.088496    31.667661 0.236
## [15,]     -4.5956090 10.565888   -30.641680     1.299086 0.316
## [16,]     -3.7256016 10.659967   -32.794314     5.974629 0.308
## [17,]    -12.7496421 17.777733   -59.090938     0.000000 0.526
## [18,]     -2.1534587  7.856944   -25.086513     3.481331 0.248
## [19,]     -2.0881613  8.620420   -24.424391     3.983888 0.222
plot(fit.npb.18$beta[,1], type = "l")

plot(fit.npb.18$beta[,2], type = "l")

plot(fit.npb.18$beta[,13], type = "l")

4.1.5.4 Try making a.phi1 = 100 and sig2inv.mu1 = 1

priors.npb.19 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 1)

fit.npb.19 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.19, interact = F)
npb.sum.19 <- summary(fit.npb.19)
npb.sum.19$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -1.93945571  7.436919   -24.014768   11.1380086 0.290
##  [2,]   -29.35356001 20.315233   -69.971480    0.0000000 0.886
##  [3,]    -0.92912783  8.390025   -23.927612   17.8853035 0.290
##  [4,]    -4.24255171 11.019911   -37.013982   11.3654517 0.352
##  [5,]    -0.02346836  7.140468   -17.342551   19.0096949 0.248
##  [6,]     0.25407353  7.350061   -14.627499   19.5678783 0.240
##  [7,]     1.57211236  8.970933   -14.205714   27.7716715 0.280
##  [8,]     5.84507537 14.003653    -9.061805   47.5658242 0.364
##  [9,]     3.04650766 10.489651   -10.050578   37.3044036 0.302
## [10,]    -0.62670622  8.574075   -20.521278   20.0444767 0.290
## [11,]    -6.73510501 13.217648   -44.885409    7.8525053 0.430
## [12,]    -2.04120535  9.872410   -29.555056   17.1520514 0.326
## [13,]    -6.11778167 12.998699   -40.797392    7.4243984 0.386
## [14,]     0.93850331  8.248327   -12.573063   24.4633603 0.286
## [15,]    -4.23998411 10.449073   -33.210287    9.0480373 0.374
## [16,]    -3.66326995 10.544825   -34.165383   10.0206015 0.338
## [17,]   -14.96004205 18.360076   -58.461572    0.4857887 0.612
## [18,]    -1.47046601 10.311921   -26.579843   20.5337795 0.336
## [19,]    -2.68239003 10.918269   -33.565613   17.0679189 0.334
plot(fit.npb.19$beta[,1], type = "l")

plot(fit.npb.19$beta[,2], type = "l")

plot(fit.npb.19$beta[,13], type = "l")

4.1.5.5 Try making a.phi1 = 100 and sig2inv.mu1 = 10

priors.npb.20 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 10)

fit.npb.20 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.20, interact = F)
npb.sum.20 <- summary(fit.npb.20)
npb.sum.20$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -1.0957585  7.253974   -19.996235    12.133194 0.260
##  [2,]    -30.7165080 20.599113   -70.129391     0.000000 0.876
##  [3,]     -1.1419660  8.278356   -23.608839    17.259197 0.272
##  [4,]     -2.8683240  7.995775   -28.839509     4.979373 0.284
##  [5,]     -0.2349687  7.355764   -15.296693    19.075479 0.248
##  [6,]      0.4651001  7.262474   -14.319838    21.345014 0.268
##  [7,]      2.3831534  9.479152   -11.090020    29.683408 0.270
##  [8,]      6.7562698 14.483801    -8.564801    48.282371 0.388
##  [9,]      1.3374750  8.602574   -11.469325    30.417517 0.262
## [10,]     -0.4535202  7.685077   -20.453290    19.872858 0.264
## [11,]     -6.8788605 13.754454   -45.639225     5.972641 0.412
## [12,]     -1.2988578  8.958916   -23.994833    13.315457 0.296
## [13,]     -6.4462750 13.404348   -43.699827     8.962873 0.420
## [14,]      0.4766060  7.309769   -14.083425    20.419289 0.254
## [15,]     -4.1468831 10.703679   -34.713864    10.198616 0.374
## [16,]     -3.7334116 11.038755   -35.748628    12.120982 0.338
## [17,]    -14.5536692 19.178035   -59.659277     0.000000 0.556
## [18,]     -1.4283640  8.830606   -24.765649    16.967138 0.288
## [19,]     -2.6752086 11.080455   -39.129165    16.421379 0.330
plot(fit.npb.20$beta[,1], type = "l")

plot(fit.npb.20$beta[,2], type = "l")

plot(fit.npb.20$beta[,13], type = "l")

4.1.5.6 Try making a.phi1 = 100 and sig2inv.mu1 = 100

priors.npb.21 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 100)

fit.npb.21 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.21, interact = F)
npb.sum.21 <- summary(fit.npb.21)
npb.sum.21$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -1.598852431  8.120628   -22.125745   14.6291093 0.350
##  [2,]  -28.636717497 20.706557   -70.367184    0.0000000 0.894
##  [3,]   -0.615990793  8.153533   -18.892649   16.9367377 0.322
##  [4,]   -4.058180165  9.708837   -33.927978   10.4658336 0.426
##  [5,]   -0.003474108  7.731554   -15.980152   23.1872641 0.322
##  [6,]    0.551072649  8.516808   -16.233849   23.5174997 0.296
##  [7,]    1.834206486  9.568153   -14.109403   29.9184057 0.338
##  [8,]    6.766290708 14.955496    -8.186632   47.3844319 0.414
##  [9,]    1.751247526  8.799965   -10.038815   26.8865154 0.290
## [10,]   -1.519519066  7.488465   -21.309451   14.6719829 0.318
## [11,]   -5.018906137 11.662643   -37.760597   11.9308543 0.434
## [12,]   -1.677747508 11.211820   -27.835487   21.5543879 0.394
## [13,]   -4.838752463 11.193045   -36.043497   10.4594162 0.430
## [14,]    0.863907073  8.691483   -13.560335   26.9581192 0.328
## [15,]   -4.327300531  9.899823   -31.088168    5.3204325 0.386
## [16,]   -3.276965016 10.419959   -31.561708   16.7571722 0.424
## [17,]  -14.099124604 18.520190   -62.348545    0.1831654 0.598
## [18,]   -0.759890932  8.582973   -18.788232   20.9963292 0.334
## [19,]   -3.076637438  9.461079   -31.613003   12.5061512 0.360
plot(fit.npb.21$beta[,1], type = "l")

plot(fit.npb.21$beta[,2], type = "l")

plot(fit.npb.21$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.21a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 100)

fit.npb.21a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = F,
                 priors = priors.npb.21a, interact = F)
npb.sum.21a <- summary(fit.npb.21a)
npb.sum.21a$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   0.0051412432 0.08469866   -0.1828181    0.2085003 0.502
##  [2,]   0.0037644320 0.09304471   -0.2130544    0.2124937 0.540
##  [3,]   0.0026560704 0.08663993   -0.1781822    0.2117353 0.530
##  [4,]   0.0013612505 0.07993577   -0.1711556    0.1806295 0.502
##  [5,]  -0.0023618224 0.08429401   -0.2003212    0.1807728 0.550
##  [6,]   0.0028257961 0.07825187   -0.1775277    0.2066858 0.486
##  [7,]  -0.0025310235 0.08453026   -0.1876556    0.1875710 0.538
##  [8,]   0.0013732366 0.08632333   -0.1878732    0.1925771 0.518
##  [9,]   0.0033625563 0.08540797   -0.1894751    0.1967161 0.520
## [10,]  -0.0049627231 0.08504732   -0.2075164    0.1905664 0.556
## [11,]  -0.0038248546 0.08145215   -0.1857765    0.1905289 0.514
## [12,]  -0.0073032570 0.08269624   -0.1918137    0.1881970 0.510
## [13,]  -0.0037368492 0.08167179   -0.1870376    0.1790697 0.504
## [14,]   0.0046757747 0.08909916   -0.1864469    0.2050904 0.574
## [15,]  -0.0021973315 0.08784817   -0.1918137    0.1872830 0.512
## [16,]   0.0019503104 0.08383080   -0.1736276    0.2080948 0.530
## [17,]   0.0077824576 0.08742294   -0.1830547    0.2071179 0.540
## [18,]   0.0001896831 0.08418336   -0.1850637    0.1905664 0.524
## [19,]   0.0063805729 0.08233081   -0.1782442    0.2015107 0.492
plot(fit.npb.21a$beta[,1], type = "l")

plot(fit.npb.21a$beta[,2], type = "l")

plot(fit.npb.21a$beta[,13], type = "l")

4.1.5.7 Try making a.phi1 = 10e5 and sig2inv.mu1 = 10e5

priors.npb.22 <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10e5, sig2inv.mu1 = 10e5)

fit.npb.22 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.22, interact = F)
npb.sum.22 <- summary(fit.npb.22)
npb.sum.22$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.035306284 0.4087304   -1.0444576    0.8734612 0.520
##  [2,]   -0.047707128 0.4067405   -0.9722595    0.8787816 0.502
##  [3,]   -0.035577082 0.4012861   -1.0507546    0.8651339 0.506
##  [4,]   -0.057457586 0.4065833   -1.0596305    0.7537023 0.526
##  [5,]   -0.044657655 0.4247126   -1.1437954    0.9025106 0.494
##  [6,]   -0.047098218 0.4545910   -1.1189796    1.0003863 0.506
##  [7,]   -0.002177765 0.4795466   -1.0556023    1.0804129 0.558
##  [8,]   -0.011803902 0.3942380   -0.9792511    0.9517408 0.498
##  [9,]   -0.033054139 0.4535246   -1.1286721    0.9043937 0.548
## [10,]   -0.043970919 0.4285394   -1.1832949    0.9432537 0.510
## [11,]   -0.034836820 0.4687112   -1.0535954    1.1372874 0.526
## [12,]   -0.055681433 0.4141023   -1.1090113    0.8180197 0.478
## [13,]   -0.061520085 0.4261508   -1.1048181    0.8484011 0.506
## [14,]   -0.014160757 0.4547002   -1.0934082    0.9126038 0.534
## [15,]   -0.047145647 0.4326528   -1.1168833    0.8238743 0.478
## [16,]   -0.048599028 0.4232817   -1.0691288    0.7790946 0.528
## [17,]   -0.070057692 0.4356459   -1.1557702    0.7609288 0.536
## [18,]   -0.020791184 0.4355400   -1.1591367    0.9563392 0.476
## [19,]   -0.040411595 0.4550934   -1.0530173    0.9865404 0.506
plot(fit.npb.22$beta[,1], type = "l")

plot(fit.npb.22$beta[,2], type = "l")

plot(fit.npb.22$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.22a <- list(alpha.pi = 8, beta.pi = 8, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10e5, sig2inv.mu1 = 10e5)

fit.npb.22a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.22a, interact = F)
npb.sum.22a <- summary(fit.npb.22a)
npb.sum.22a$main.effects
##        Posterior Mean           SD 95% CI Lower 95% CI Upper   PIP
##  [1,]  0.000006182868 0.0008619421 -0.001979458  0.002071416 0.520
##  [2,] -0.000052064014 0.0008040348 -0.001874284  0.001670370 0.524
##  [3,] -0.000046371401 0.0008962762 -0.002141472  0.001883780 0.522
##  [4,]  0.000039298564 0.0009298157 -0.002246873  0.001975836 0.522
##  [5,] -0.000017857840 0.0008660585 -0.002082111  0.001861066 0.546
##  [6,] -0.000076962013 0.0009002153 -0.002196562  0.001778247 0.512
##  [7,] -0.000017274586 0.0008623192 -0.001970859  0.002013509 0.528
##  [8,] -0.000042419780 0.0008878686 -0.001940788  0.001864174 0.528
##  [9,] -0.000037825550 0.0008999260 -0.002129650  0.001889525 0.542
## [10,] -0.000071637744 0.0008846355 -0.002220330  0.001768216 0.506
## [11,]  0.000018945719 0.0009117751 -0.002224579  0.002037273 0.520
## [12,]  0.000076549490 0.0007989603 -0.002009522  0.001835409 0.470
## [13,]  0.000065717471 0.0008749839 -0.001984107  0.001916398 0.524
## [14,] -0.000030690991 0.0008669480 -0.002107817  0.001895524 0.550
## [15,]  0.000019370894 0.0009155382 -0.001984431  0.001988277 0.542
## [16,] -0.000036712651 0.0009179785 -0.002146234  0.001986838 0.502
## [17,] -0.000028941612 0.0008525617 -0.001949226  0.001762325 0.512
## [18,]  0.000010433393 0.0009188860 -0.002237139  0.002022782 0.514
## [19,] -0.000023885689 0.0008964864 -0.002187121  0.001901227 0.540
plot(fit.npb.22a$beta[,1], type = "l")

plot(fit.npb.22a$beta[,2], type = "l")

plot(fit.npb.22a$beta[,13], type = "l")

4.1.6 Set alpha.pi and beta.pi to 5, readjust a.phi1 and sig2inv.mu1

Set alpha.pi and beta.pi to 5, rather than 8, and try adjusting a.phi1 and sig2inv.mu1

4.1.6.1 Try making a.phi1 = 10 and sig2inv.mu1 = 1

priors.npb.23 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 1)

fit.npb.23 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.23, interact = F)
npb.sum.23 <- summary(fit.npb.23)
npb.sum.23$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -2.2669275  6.818737   -23.985602     7.226331 0.284
##  [2,]    -23.8968298 21.256605   -70.757009     0.000000 0.802
##  [3,]     -1.3764332  5.783617   -18.251351     6.223349 0.220
##  [4,]     -3.1824556  8.462735   -27.548227     3.788901 0.308
##  [5,]     -1.1498466  5.614243   -17.671669     6.601940 0.232
##  [6,]     -1.0777643  6.557314   -19.429210    11.728310 0.250
##  [7,]      0.3905512  8.642639   -14.259111    24.121953 0.232
##  [8,]      3.6174570 12.688717    -9.551844    45.028736 0.256
##  [9,]      0.8883675  8.857065   -13.113582    23.017336 0.236
## [10,]     -0.5558551  5.804998   -14.312065     8.953615 0.212
## [11,]     -4.8874443 10.554609   -35.472399     5.931729 0.372
## [12,]     -2.0602220  7.428268   -24.099109     8.386962 0.280
## [13,]     -5.3142433 10.814921   -40.185742     1.913820 0.358
## [14,]     -0.6998181  6.073537   -16.303471    11.832843 0.236
## [15,]     -2.9213200  7.868612   -26.124243     4.919153 0.308
## [16,]     -3.8903774 10.524390   -35.856330     8.674976 0.354
## [17,]    -11.2928049 17.090541   -56.513393     1.655472 0.514
## [18,]     -1.8947305  7.612136   -22.642968     8.646449 0.276
## [19,]     -3.2289844  8.910317   -28.874166     3.991413 0.294
plot(fit.npb.23$beta[,1], type = "l")

plot(fit.npb.23$beta[,2], type = "l")

plot(fit.npb.23$beta[,13], type = "l")

Same priors, but now with scaleY = F

priors.npb.23a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 1)

fit.npb.23a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.23a, interact = F)
npb.sum.23a <- summary(fit.npb.23a)
npb.sum.23a$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]  -0.0007196103 0.2964622   -0.6553441    0.7020867 0.510
##  [2,]  -0.0120527174 0.2632085   -0.6428513    0.5712860 0.498
##  [3,]  -0.0003302213 0.2964493   -0.6474469    0.6524802 0.538
##  [4,]   0.0038542799 0.2704418   -0.6072693    0.5571160 0.502
##  [5,]   0.0073328019 0.2603254   -0.5727369    0.5508989 0.472
##  [6,]   0.0149783637 0.2940539   -0.6596654    0.7206298 0.482
##  [7,]  -0.0148465978 0.3064579   -0.6605561    0.6640202 0.534
##  [8,]   0.0042084310 0.3055809   -0.6986858    0.6806700 0.502
##  [9,]   0.0040830043 0.2701805   -0.5430343    0.6331375 0.480
## [10,]   0.0070490226 0.2993666   -0.6461819    0.7809095 0.516
## [11,]   0.0104681999 0.2989888   -0.6058662    0.7507703 0.524
## [12,]  -0.0009639157 0.2873621   -0.6392656    0.6505214 0.512
## [13,]  -0.0202041392 0.2846702   -0.6886220    0.5432586 0.532
## [14,]   0.0097133647 0.2979735   -0.6236524    0.7223282 0.506
## [15,]   0.0268837151 0.3114360   -0.6323325    0.7621870 0.508
## [16,]  -0.0019335793 0.2594860   -0.6081490    0.6522747 0.462
## [17,]  -0.0144824164 0.2984423   -0.7850978    0.6364156 0.514
## [18,]   0.0088153563 0.2645159   -0.5405364    0.7058493 0.468
## [19,]   0.0197624609 0.3104259   -0.6356867    0.7961824 0.490
plot(fit.npb.23a$beta[,1], type = "l")

plot(fit.npb.23a$beta[,2], type = "l")

plot(fit.npb.23a$beta[,13], type = "l")

4.1.6.2 Try making a.phi1 = 10 and sig2inv.mu1 = 10

priors.npb.24 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 10)

fit.npb.24 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.24, interact = F)
npb.sum.24 <- summary(fit.npb.24)
npb.sum.24$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -1.1744090  6.356825   -21.365240    6.5471919 0.180
##  [2,]    -25.9222430 21.582287   -70.704025    0.0000000 0.764
##  [3,]     -0.8719158  5.355631   -17.380790    5.0388036 0.154
##  [4,]     -2.6852993  8.750953   -31.977107    0.8560601 0.214
##  [5,]     -0.3376489  4.707140   -13.656701    4.3798985 0.122
##  [6,]     -0.1516117  6.011091   -12.307358   11.9384392 0.158
##  [7,]      0.2033702  5.280983    -9.111321   13.6905185 0.150
##  [8,]      4.1573659 12.676688    -5.377147   45.2223780 0.212
##  [9,]      0.8289965  7.574134   -10.641479   24.7475649 0.164
## [10,]     -0.6102648  6.032846   -16.157754   10.9900904 0.156
## [11,]     -4.5080221 10.836041   -37.483169    2.7381582 0.274
## [12,]     -1.5336841  7.601180   -23.358433    3.3745605 0.192
## [13,]     -4.5377605 11.415565   -39.490858    0.9601621 0.280
## [14,]     -0.1297684  6.641396   -12.858573   18.9852703 0.178
## [15,]     -3.8896588 10.341413   -37.505833    2.9042601 0.282
## [16,]     -3.6972612  9.804580   -36.883240    1.6065287 0.230
## [17,]    -11.4272777 18.220939   -58.972882    0.0000000 0.434
## [18,]     -1.5892077  7.539632   -22.211002    8.4387392 0.204
## [19,]     -2.7649340  9.661202   -36.632517    3.6407285 0.210
plot(fit.npb.24$beta[,1], type = "l")

plot(fit.npb.24$beta[,2], type = "l")

plot(fit.npb.24$beta[,13], type = "l")

4.1.6.3 Try making a.phi1 = 10 and sig2inv.mu1 = 100

priors.npb.25 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10, sig2inv.mu1 = 100)

fit.npb.25 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.25, interact = F)
npb.sum.25 <- summary(fit.npb.25)
npb.sum.25$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -1.66648315  7.226089   -26.871111    6.1849869 0.200
##  [2,]   -28.80748355 20.864602   -67.769965    0.0000000 0.838
##  [3,]    -1.42558012  6.522784   -22.771262    4.4281229 0.192
##  [4,]    -2.76921892  8.454176   -31.199743    1.5756317 0.216
##  [5,]    -0.28525201  4.552973   -10.345389    7.3140781 0.144
##  [6,]    -0.02273501  5.207522   -10.280686   12.2972152 0.128
##  [7,]     0.42838036  6.509331    -8.832926   15.5831741 0.146
##  [8,]     3.61899464 11.611365    -4.922887   41.1492615 0.204
##  [9,]     0.65651537  6.752541    -9.044913   17.4790875 0.154
## [10,]    -0.28000249  4.486445   -11.157330    7.2672729 0.164
## [11,]    -5.83525887 12.531493   -42.110405    2.3667362 0.306
## [12,]    -2.20787228  8.966495   -29.191497    3.2932711 0.220
## [13,]    -7.12779870 14.099981   -48.626023    1.9446005 0.344
## [14,]     0.36657831  8.642415   -15.834819   24.3479811 0.194
## [15,]    -4.05853291 10.317329   -35.917289    0.7062024 0.242
## [16,]    -3.19567962  9.827098   -34.041145    4.7932942 0.244
## [17,]   -12.25439359 18.517567   -58.794296    0.0000000 0.434
## [18,]    -1.90909272  8.525734   -27.250940    5.5485470 0.204
## [19,]    -3.46970786 10.432867   -37.223133    3.3684512 0.248
plot(fit.npb.25$beta[,1], type = "l")

plot(fit.npb.25$beta[,2], type = "l")

plot(fit.npb.25$beta[,13], type = "l")

4.1.6.4 Try making a.phi1 = 100 and sig2inv.mu1 = 1

priors.npb.26 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 1)

fit.npb.26 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.26, interact = F)
npb.sum.26 <- summary(fit.npb.26)
npb.sum.26$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -1.2621545  7.405051   -22.005981    11.964441 0.262
##  [2,]    -31.8730363 20.471567   -69.939915     0.000000 0.862
##  [3,]     -0.7504150  8.021987   -20.950051    16.075527 0.282
##  [4,]     -4.1027594 11.475524   -41.354956     9.531501 0.336
##  [5,]      0.1652421  5.658890   -14.033363    14.483537 0.236
##  [6,]      0.7633474  7.314781   -14.251274    18.801390 0.296
##  [7,]      2.1124896  9.166049   -11.404178    28.577810 0.266
##  [8,]      7.6878235 15.418804    -7.034124    50.178836 0.410
##  [9,]      2.7160023  9.653674    -9.890912    31.315221 0.314
## [10,]      0.3323067  8.778377   -17.837627    18.701679 0.288
## [11,]     -6.5489320 13.681433   -44.351904     6.346462 0.378
## [12,]     -2.2381381  9.969420   -30.264782    15.145125 0.320
## [13,]     -4.4238949 11.719498   -37.106223     9.696419 0.354
## [14,]      1.1390514  7.813899   -13.346768    22.660832 0.264
## [15,]     -3.6777909 10.307871   -34.607732     9.943842 0.330
## [16,]     -1.4154880 10.777865   -31.506473    19.824269 0.322
## [17,]    -16.4797330 20.840367   -63.168555     1.184975 0.564
## [18,]     -0.6581118  8.235614   -23.385310    18.644877 0.318
## [19,]     -1.6463141  9.357480   -30.512594    13.744361 0.288
plot(fit.npb.26$beta[,1], type = "l")

plot(fit.npb.26$beta[,2], type = "l")

plot(fit.npb.26$beta[,13], type = "l")

4.1.6.5 Try making a.phi1 = 100 and sig2inv.mu1 = 10

priors.npb.27 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 10)

fit.npb.27 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.27, interact = F)
npb.sum.27 <- summary(fit.npb.27)
npb.sum.27$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -1.01163147  6.974847   -23.305300    12.063268 0.248
##  [2,]   -30.66127687 20.760120   -70.934734     0.000000 0.872
##  [3,]    -1.00436024  7.982394   -24.827645    15.695069 0.270
##  [4,]    -3.22484204 10.981124   -35.503109    10.899975 0.304
##  [5,]     0.65483898  7.215895   -14.305391    21.422911 0.226
##  [6,]     1.05202716  8.610743   -13.386552    26.288227 0.240
##  [7,]     1.04745993  9.161936   -18.039125    24.969603 0.270
##  [8,]     7.61348945 15.349643    -9.331034    48.223341 0.388
##  [9,]     2.22017934  9.675796   -14.805064    31.309280 0.256
## [10,]     0.06297748  7.844041   -17.031138    18.704314 0.260
## [11,]    -6.49365289 13.166477   -44.745196     5.150326 0.368
## [12,]    -1.23067543 10.029990   -29.537828    21.137705 0.272
## [13,]    -5.58891305 13.398458   -44.759448     6.707673 0.390
## [14,]     0.48784423  9.658446   -21.374144    27.318140 0.246
## [15,]    -3.70084647  9.962155   -34.880874     3.908209 0.288
## [16,]    -1.89826339 10.065355   -26.207320    15.512113 0.304
## [17,]   -14.86352388 19.333278   -60.535336     0.000000 0.540
## [18,]    -1.82229314  8.936614   -29.039057    12.853870 0.270
## [19,]    -3.11206478 11.992936   -39.216810    15.476387 0.326
plot(fit.npb.27$beta[,1], type = "l")

plot(fit.npb.27$beta[,2], type = "l")

plot(fit.npb.27$beta[,13], type = "l")

4.1.6.6 Try making a.phi1 = 100 and sig2inv.mu1 = 100

priors.npb.28 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 100)

fit.npb.28 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.28, interact = F)
npb.sum.28 <- summary(fit.npb.28)
npb.sum.28$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -1.5967933  7.684725   -23.845516   12.9857927 0.262
##  [2,]    -30.3509594 21.129884   -69.873352    0.0000000 0.848
##  [3,]     -1.1081974  7.846097   -22.418736   14.0079554 0.260
##  [4,]     -3.1540371  9.566356   -32.632348    8.9516887 0.292
##  [5,]     -0.1257847  6.959253   -16.825701   18.5528850 0.254
##  [6,]      0.5449194  7.866129   -13.596188   22.6291587 0.288
##  [7,]      1.1863816  9.553280   -13.924976   27.4139131 0.288
##  [8,]      8.7649384 17.034173   -10.170230   53.1274332 0.444
##  [9,]      3.5466042 11.407061    -9.182936   41.9792743 0.328
## [10,]     -0.3097942  7.638762   -19.690696   17.8438100 0.274
## [11,]     -5.9581373 12.160559   -40.851977    6.9929063 0.382
## [12,]     -1.6259110 10.051129   -26.222054   19.4293422 0.316
## [13,]     -6.1755473 13.010573   -42.585762    7.8574922 0.394
## [14,]      1.7462565 10.147823   -17.436863   33.7645626 0.272
## [15,]     -3.6360794 10.305510   -36.022150    9.1708419 0.324
## [16,]     -1.4794849 11.214146   -30.263661   21.3994828 0.330
## [17,]    -17.2136283 20.264874   -63.256164    0.4187659 0.618
## [18,]     -1.0341528  9.680756   -27.793811   19.1325054 0.306
## [19,]     -2.1742274 11.416120   -33.177502   19.8048518 0.328
plot(fit.npb.28$beta[,1], type = "l")

plot(fit.npb.28$beta[,2], type = "l")

plot(fit.npb.28$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.28a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 100, sig2inv.mu1 = 100)

fit.npb.28a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = F,
                 priors = priors.npb.28a, interact = F)
npb.sum.28a <- summary(fit.npb.28a)
npb.sum.28a$main.effects
##       Posterior Mean         SD 95% CI Lower 95% CI Upper   PIP
##  [1,]  0.00304819944 0.07585417   -0.1791721    0.1588592 0.504
##  [2,]  0.00368920822 0.08356244   -0.1832763    0.1953595 0.484
##  [3,]  0.00230653550 0.08395801   -0.1981486    0.1913326 0.518
##  [4,]  0.00349433034 0.08395506   -0.1834283    0.2139654 0.484
##  [5,] -0.00392341030 0.07385651   -0.1812940    0.1651630 0.480
##  [6,] -0.00008981685 0.08121901   -0.1920336    0.1843800 0.484
##  [7,]  0.00700464294 0.08178275   -0.1759053    0.2133049 0.480
##  [8,]  0.00157109764 0.07703858   -0.1735443    0.1754262 0.476
##  [9,]  0.01105053383 0.07811738   -0.1558621    0.2026439 0.496
## [10,]  0.00674268321 0.08007177   -0.1618587    0.1982301 0.484
## [11,]  0.00070700654 0.07406451   -0.1659681    0.1679017 0.482
## [12,]  0.00148550141 0.08263396   -0.1969318    0.1810100 0.502
## [13,]  0.00514617985 0.08217526   -0.1626479    0.1919755 0.514
## [14,] -0.00020879633 0.08166007   -0.1723812    0.1923978 0.502
## [15,] -0.00140959112 0.07936366   -0.1735443    0.1720113 0.488
## [16,]  0.00212519609 0.07457426   -0.1619751    0.1777630 0.494
## [17,]  0.00486790268 0.07542233   -0.1648815    0.1680491 0.516
## [18,] -0.00225180822 0.08085503   -0.1661710    0.1791710 0.516
## [19,]  0.00332850383 0.07821127   -0.1611759    0.1848179 0.520
plot(fit.npb.28a$beta[,1], type = "l")

plot(fit.npb.28a$beta[,2], type = "l")

plot(fit.npb.28a$beta[,13], type = "l")

4.1.6.7 Try making a.phi1 = 10e5 and sig2inv.mu1 = 10e5

priors.npb.29 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10e5, sig2inv.mu1 = 10e5)

fit.npb.29 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.29, interact = F)
npb.sum.29 <- summary(fit.npb.29)
npb.sum.29$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]   -0.006677226 0.4192719   -0.9510673    0.9945879 0.492
##  [2,]   -0.069118614 0.4012581   -1.1079014    0.7428316 0.484
##  [3,]   -0.038108583 0.4021976   -1.0728872    0.8241823 0.506
##  [4,]   -0.048976955 0.4261434   -1.0467512    0.8276469 0.532
##  [5,]   -0.053049503 0.4064875   -1.0925019    0.7353015 0.498
##  [6,]   -0.018892737 0.4176082   -0.9341080    0.9094812 0.532
##  [7,]   -0.055467763 0.4250969   -1.0942571    0.8192562 0.490
##  [8,]   -0.007563029 0.4117104   -0.9800507    0.8741712 0.548
##  [9,]   -0.032509036 0.4138104   -0.9341025    0.9188171 0.498
## [10,]   -0.050876901 0.4142152   -1.0858341    0.8606469 0.548
## [11,]   -0.043178598 0.4017999   -1.0177024    0.8189672 0.540
## [12,]   -0.049915027 0.4154000   -0.9444616    0.8452286 0.518
## [13,]   -0.044566791 0.4348905   -0.9736428    0.9201996 0.542
## [14,]   -0.053132877 0.4085925   -0.9758844    0.8820159 0.482
## [15,]   -0.046278234 0.3940232   -0.9180320    0.8092658 0.500
## [16,]   -0.049514433 0.4578257   -1.1040893    0.9481966 0.534
## [17,]   -0.070271101 0.4259853   -1.0340079    0.8111130 0.532
## [18,]   -0.042616463 0.3950922   -1.0527976    0.8091935 0.482
## [19,]   -0.039944456 0.4197640   -1.0431318    0.8395898 0.502
plot(fit.npb.29$beta[,1], type = "l")

plot(fit.npb.29$beta[,2], type = "l")

plot(fit.npb.29$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.29a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 10e5, sig2inv.mu1 = 10e5)

fit.npb.29a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.29a, interact = F)
npb.sum.29a <- summary(fit.npb.29a)
npb.sum.29a$main.effects
##         Posterior Mean           SD 95% CI Lower 95% CI Upper   PIP
##  [1,] -0.0000029929547 0.0008091636 -0.001938430  0.001936734 0.522
##  [2,]  0.0000348936641 0.0007978201 -0.001778601  0.001874876 0.534
##  [3,]  0.0000050700913 0.0008528772 -0.001883928  0.001829383 0.502
##  [4,] -0.0000006592784 0.0008320206 -0.001783612  0.002043350 0.496
##  [5,]  0.0000070091150 0.0007874591 -0.001657262  0.001945618 0.496
##  [6,] -0.0000142948030 0.0008108816 -0.001693869  0.001954015 0.512
##  [7,]  0.0000299715926 0.0008253950 -0.001715913  0.002258297 0.506
##  [8,]  0.0000180048470 0.0008612665 -0.002026700  0.001880374 0.524
##  [9,]  0.0000777242366 0.0007818032 -0.001610301  0.002026431 0.480
## [10,] -0.0000037155238 0.0008344709 -0.001948587  0.001890875 0.544
## [11,]  0.0000539999709 0.0008191653 -0.001793009  0.001999934 0.506
## [12,] -0.0000143277900 0.0007643464 -0.001657805  0.001806774 0.476
## [13,]  0.0000019952710 0.0008151975 -0.001850042  0.001889987 0.486
## [14,] -0.0000099848829 0.0008049083 -0.001912550  0.001996068 0.496
## [15,]  0.0000398545360 0.0007695658 -0.001654104  0.001944979 0.506
## [16,]  0.0000044317518 0.0008533496 -0.001914691  0.001868051 0.542
## [17,]  0.0000318213211 0.0007992530 -0.001655477  0.001985471 0.492
## [18,]  0.0000436334011 0.0007482470 -0.001591704  0.001803765 0.468
## [19,]  0.0000257047232 0.0007601775 -0.001762099  0.001869755 0.470
plot(fit.npb.29a$beta[,1], type = "l")

plot(fit.npb.29a$beta[,2], type = "l")

plot(fit.npb.29a$beta[,13], type = "l")

4.1.7 Set alpha.pi and beta.pi to 5, set a.phi1 = 1, and adjust sig2inv.mu1

4.1.7.1 Try making a.phi1 = 1 and sig2inv.mu1 = 10

priors.npb.30 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 10)

fit.npb.30 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.30, interact = F)
npb.sum.30 <- summary(fit.npb.30)
npb.sum.30$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]     -1.7228556  6.185898   -20.101281    1.8610950 0.200
##  [2,]    -19.0050772 20.190536   -65.928054    0.0000000 0.684
##  [3,]     -1.3136830  5.208024   -19.738044    1.6554852 0.172
##  [4,]     -3.1161119  8.173765   -28.954626    0.2456618 0.236
##  [5,]     -0.4698456  4.195973   -11.667282    3.9323218 0.144
##  [6,]     -0.4558571  5.069872   -12.863175    6.0642046 0.160
##  [7,]      0.3000572  6.026863    -8.381134   17.8194398 0.154
##  [8,]      2.1301516 10.686663    -8.617136   41.4826373 0.194
##  [9,]      0.2229029  6.518747    -9.693930   10.3608104 0.144
## [10,]     -1.3429782  5.071480   -16.423261    1.5875183 0.178
## [11,]     -3.9985155  9.784726   -35.178073    2.2364369 0.278
## [12,]     -2.7480595  7.784423   -27.773997    1.6691937 0.242
## [13,]     -4.3861962  9.853896   -35.243232    0.0000000 0.288
## [14,]     -0.6861690  6.209898   -14.500452    2.4538156 0.170
## [15,]     -3.5042988  8.989247   -32.428574    0.2456618 0.258
## [16,]     -3.9430599  9.088170   -32.402623    0.4364748 0.284
## [17,]     -9.4313789 16.160430   -50.994577    0.4967513 0.420
## [18,]     -2.0118603  6.588875   -19.891292    0.0000000 0.212
## [19,]     -1.6924124  6.275389   -19.649915    0.8215062 0.170
plot(fit.npb.30$beta[,1], type = "l")

plot(fit.npb.30$beta[,2], type = "l")

plot(fit.npb.30$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.30a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 10)

fit.npb.30a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.30a, interact = F)
npb.sum.30a <- summary(fit.npb.30a)
npb.sum.30a$main.effects
##       Posterior Mean       SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    0.074119118 1.363934    -2.308498     2.756957 0.468
##  [2,]    0.010494643 1.326295    -2.420825     2.808571 0.512
##  [3,]    0.063236533 1.262993    -2.227545     2.588655 0.482
##  [4,]    0.004503723 1.485491    -2.346604     2.475782 0.488
##  [5,]   -0.018773342 1.337000    -2.403715     2.326143 0.510
##  [6,]   -0.014432961 1.298491    -2.198557     2.679880 0.476
##  [7,]    0.098224373 1.286374    -1.921173     2.803585 0.476
##  [8,]   -0.046045737 1.194057    -2.256074     2.550880 0.502
##  [9,]    0.060598256 1.322770    -2.188775     2.812985 0.516
## [10,]    0.103873445 1.319422    -2.204537     2.905311 0.488
## [11,]    0.067135740 1.344225    -2.214774     1.920478 0.530
## [12,]    0.064489504 1.228626    -2.179235     2.507171 0.496
## [13,]   -0.038101943 1.473007    -3.643061     2.745544 0.480
## [14,]    0.056993085 1.264079    -2.323791     3.011037 0.496
## [15,]    0.043757160 1.258120    -2.222910     2.656063 0.504
## [16,]    0.121804354 1.072164    -1.878659     2.567970 0.468
## [17,]   -0.020926185 1.240334    -3.079505     2.507171 0.492
## [18,]    0.072878661 1.423697    -2.532435     2.588655 0.480
## [19,]   -0.032569195 1.034985    -1.843805     1.854254 0.456
plot(fit.npb.30a$beta[,1], type = "l")

plot(fit.npb.30a$beta[,2], type = "l")

plot(fit.npb.30a$beta[,13], type = "l")

4.1.7.2 Try making a.phi1 = 1 and sig2inv.mu1 = 100

priors.npb.31 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 100)

fit.npb.31 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.31, interact = F)
npb.sum.31 <- summary(fit.npb.31)
npb.sum.31$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -2.30889214  6.784902   -24.771436  0.000000000 0.200
##  [2,]   -20.92802151 20.179046   -62.805035  0.000000000 0.732
##  [3,]    -1.08965652  4.397055   -12.799255  0.000000000 0.158
##  [4,]    -3.37844361  8.694089   -32.124019  0.003170527 0.252
##  [5,]    -0.93522724  4.686011   -13.243451  0.188809267 0.156
##  [6,]    -0.54554790  5.140308   -12.597600  6.688217638 0.164
##  [7,]    -0.19526708  5.578431   -11.915487 12.222610797 0.184
##  [8,]     1.75217377  9.034742    -7.926483 33.065635224 0.174
##  [9,]     0.09622656  5.085413    -8.919991 12.543905370 0.146
## [10,]    -1.07175896  6.258691   -16.323661  6.033725382 0.192
## [11,]    -4.87259495 10.790576   -36.999806  0.003170527 0.310
## [12,]    -2.65268645  7.873357   -25.620569  0.000000000 0.238
## [13,]    -4.77325818 10.868420   -37.171170  0.038797733 0.314
## [14,]    -0.70896380  5.283107   -12.090274  5.821940112 0.172
## [15,]    -3.02747698  8.213266   -28.371188  0.538732726 0.266
## [16,]    -3.11612392  8.292508   -28.674739  0.000000000 0.260
## [17,]   -11.12370533 17.011735   -56.854870  0.000000000 0.454
## [18,]    -2.24537034  7.309293   -26.808510  0.108378804 0.214
## [19,]    -2.18288162  6.462647   -24.682852  0.038797733 0.220
plot(fit.npb.31$beta[,1], type = "l")

plot(fit.npb.31$beta[,2], type = "l")

plot(fit.npb.31$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.31a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 100)

fit.npb.31a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.31a, interact = F)
npb.sum.31a <- summary(fit.npb.31a)
npb.sum.31a$main.effects
##       Posterior Mean       SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -0.11200176 1.845975    -3.594599     2.933795 0.524
##  [2,]    -0.06084745 1.446065    -3.230574     2.402961 0.544
##  [3,]    -0.03186257 1.609605    -3.502165     2.666465 0.494
##  [4,]    -0.05796699 1.565341    -3.072019     2.393704 0.542
##  [5,]    -0.03733648 1.353984    -2.690057     2.838679 0.500
##  [6,]    -0.10830473 1.893346    -5.062828     2.871385 0.532
##  [7,]    -0.09578467 1.428884    -3.525617     2.485758 0.504
##  [8,]    -0.10651158 1.903310    -3.568060     2.641535 0.518
##  [9,]    -0.07269353 1.475613    -3.204855     2.616656 0.474
## [10,]    -0.16545997 1.853086    -3.372041     2.382987 0.518
## [11,]    -0.18876057 1.724617    -3.552851     2.094670 0.478
## [12,]    -0.09518493 1.543222    -3.608179     2.813356 0.550
## [13,]    -0.08358490 1.709686    -3.763812     2.924885 0.526
## [14,]     0.01514462 1.481139    -3.398926     3.719411 0.490
## [15,]    -0.16326693 1.731134    -3.794918     2.855899 0.508
## [16,]    -0.02413426 1.722385    -3.230866     3.497753 0.518
## [17,]    -0.10137685 1.957256    -3.491128     2.648940 0.482
## [18,]     0.02946311 1.412729    -2.854176     2.989497 0.522
## [19,]    -0.02930316 1.780401    -3.267473     2.635170 0.524
plot(fit.npb.31a$beta[,1], type = "l")

plot(fit.npb.31a$beta[,2], type = "l")

plot(fit.npb.31a$beta[,13], type = "l")

4.1.7.3 Try making a.phi1 = 1 and sig2inv.mu1 = 10e5

priors.npb.32 <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 10e5)

fit.npb.32 <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = TRUE,
                 priors = priors.npb.32, interact = F)
npb.sum.32 <- summary(fit.npb.32)
npb.sum.32$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    -1.38446846  7.022531   -20.332156  10.17422656 0.226
##  [2,]   -20.62264930 20.077258   -65.089736   0.08685412 0.720
##  [3,]    -1.06435183  4.917391   -14.901425   2.53796662 0.166
##  [4,]    -3.06679682  8.528214   -27.692337   2.01332470 0.250
##  [5,]    -0.84960024  5.017413   -16.335625   5.41330050 0.146
##  [6,]    -0.38439458  5.549323   -13.393940   8.34579300 0.164
##  [7,]     0.03984075  5.288480   -12.413025  10.51676007 0.148
##  [8,]     2.37295219 10.314854    -6.581860  39.89358420 0.170
##  [9,]     0.54030411  6.489239    -9.616982  16.48935580 0.166
## [10,]    -0.64553010  5.324006   -14.722058   7.73864664 0.168
## [11,]    -5.08187604 11.058515   -36.058333   2.46512327 0.308
## [12,]    -3.02452916  8.638986   -31.762681   0.98918597 0.226
## [13,]    -4.79160611 10.602028   -37.257283   0.42522579 0.290
## [14,]    -0.47361111  5.411474   -12.858097  11.24888674 0.170
## [15,]    -3.14471481  8.448053   -30.112915   2.14570556 0.232
## [16,]    -2.96886714  9.279894   -34.730379   4.91918466 0.250
## [17,]   -10.16885803 16.029957   -51.528538   0.26003816 0.430
## [18,]    -2.44958126  8.455490   -28.831080   4.40927059 0.208
## [19,]    -2.15767429  7.679553   -24.673937   7.73864664 0.214
plot(fit.npb.32$beta[,1], type = "l")

plot(fit.npb.32$beta[,2], type = "l")

plot(fit.npb.32$beta[,13], type = "l")

Same priors, but with scaleY = F

priors.npb.32a <- list(alpha.pi = 5, beta.pi = 5, alpha.pi2 = 9, beta.pi2 = 1,
                     a.phi1 = 1, sig2inv.mu1 = 10e5)

fit.npb.32a <- npb(niter = 1000, nburn = 500, X = X.scaled, Y = Y, W = W.scaled2,
                 scaleY = FALSE,
                 priors = priors.npb.32a, interact = F)
npb.sum.32a <- summary(fit.npb.32a)
npb.sum.32a$main.effects
##       Posterior Mean       SD 95% CI Lower 95% CI Upper   PIP
##  [1,]    0.025209362 1.294861    -2.347621     2.420512 0.514
##  [2,]   -0.049016935 1.341006    -2.409473     2.263052 0.494
##  [3,]    0.028896975 1.381161    -1.960469     2.326273 0.518
##  [4,]   -0.038418735 1.219774    -2.483217     1.975696 0.518
##  [5,]   -0.023985814 1.355683    -2.392673     2.356443 0.474
##  [6,]    0.030791522 1.302866    -2.599089     2.388541 0.478
##  [7,]   -0.152837433 1.720482    -2.534500     1.739969 0.486
##  [8,]   -0.040706452 1.508916    -2.671380     2.258912 0.492
##  [9,]    0.092356015 1.371136    -2.119311     2.417871 0.480
## [10,]   -0.075665075 1.216456    -2.298597     2.012301 0.496
## [11,]    0.029002313 1.263324    -2.564436     2.791453 0.486
## [12,]   -0.099399846 1.027492    -2.148132     2.105168 0.436
## [13,]    0.152498936 1.428497    -2.158852     3.180683 0.498
## [14,]    0.062607859 1.158282    -1.939363     2.221776 0.482
## [15,]   -0.054083611 1.694672    -2.491200     2.768224 0.496
## [16,]   -0.064899159 1.629433    -3.019469     2.619183 0.514
## [17,]   -0.099680901 1.190646    -2.386633     1.873959 0.502
## [18,]   -0.008919667 1.258053    -2.299866     2.293204 0.496
## [19,]    0.092523129 1.529017    -2.054268     2.607070 0.494
plot(fit.npb.32a$beta[,1], type = "l")

plot(fit.npb.32a$beta[,2], type = "l")

plot(fit.npb.32a$beta[,13], type = "l")

4.2 Fit the NPB model

Below I’ve used the 24th set of priors and set scaleY = T

The priors are as follows: r priors.npb.24

Note that this version of the model does not include gest_age_w.

priors.npb <- priors.npb.24

fit.npb <- npb(niter = 10000, nburn = 2500, X = X.scaled, Y = Y, W = W.scaled2,
               scaleY = TRUE,
               priors = priors.npb, interact = TRUE, XWinteract = TRUE)
save(fit.npb, file = here::here("Results", "NPB_Birth_Weight.rdata"))

# load(here::here("Results", "NPB_Birth_Weight.rdata"))
npb.sum <- summary(fit.npb)

4.2.1 First, main effect regression coefficients with PIPs

npb.sum$main.effects
##       Posterior Mean        SD 95% CI Lower 95% CI Upper       PIP
##  [1,]     -1.8661265  6.914319   -22.957048     6.049375 0.2281333
##  [2,]    -24.9299099 21.166516   -69.288505     0.000000 0.7853333
##  [3,]     -1.1551286  6.033682   -18.834106     8.098385 0.2098667
##  [4,]     -3.1530965  9.092592   -32.056663     3.639703 0.2676000
##  [5,]     -0.5165156  5.370372   -14.232302     9.672763 0.1869333
##  [6,]     -0.2570868  5.721849   -14.272370    12.276842 0.1908000
##  [7,]      0.5523885  7.135228   -11.182066    19.879487 0.1908000
##  [8,]      3.8083493 12.832219    -7.874556    47.771370 0.2461333
##  [9,]      0.7159346  7.048218   -10.640137    20.509450 0.1908000
## [10,]     -0.8666888  5.859006   -17.095504     7.869133 0.1969333
## [11,]     -4.7334937 11.084603   -39.038464     2.605766 0.3028000
## [12,]     -2.3870316  8.563241   -29.073269     6.775860 0.2480000
## [13,]     -5.3488304 11.913380   -42.129727     1.969438 0.3268000
## [14,]     -0.4622358  6.634494   -15.841082    12.762775 0.2030667
## [15,]     -3.3172436  8.902109   -31.384447     3.383887 0.2705333
## [16,]     -3.0570276  9.676074   -32.593852     7.114270 0.2761333
## [17,]    -12.3144543 18.075048   -57.590304     0.000000 0.4890667
## [18,]     -1.7492751  7.631214   -23.981519     7.848185 0.2302667
## [19,]     -2.4616124  9.410943   -31.212532     7.106345 0.2552000
selected_exp <- which(npb.sum$main.effects[,5] >= 0.5)
selected_exp
## [1] 2
#' Which variables are these?
exp_names <- colnames(X)[selected_exp]
exp_names
## [1] "mean_o3"

4.2.2 Interactions

Next, all of the interactions between exposures or between exposures and covariates

npb.sum$interactions
##        Posterior Mean         SD 95% CI Lower 95% CI Upper          PIP
##   [1,]  0.00207935185  0.3922357            0            0 0.0026666667
##   [2,] -0.01092256325  0.8419317            0            0 0.0030666667
##   [3,]  0.00406659348  0.2943407            0            0 0.0012000000
##   [4,]  0.02102702601  0.8433957            0            0 0.0018666667
##   [5,] -0.00120801985  0.2514967            0            0 0.0012000000
##   [6,] -0.00434257270  0.5189284            0            0 0.0025333333
##   [7,]  0.00162709402  0.2961235            0            0 0.0013333333
##   [8,] -0.00828639526  0.7148152            0            0 0.0025333333
##   [9,] -0.04316211873  1.1765333            0            0 0.0040000000
##  [10,] -0.01240167260  0.4834403            0            0 0.0017333333
##  [11,] -0.00816117136  0.6874690            0            0 0.0024000000
##  [12,] -0.00174916587  0.3175548            0            0 0.0020000000
##  [13,] -0.03365644311  0.9396474            0            0 0.0025333333
##  [14,] -0.00081702879  0.3060190            0            0 0.0014666667
##  [15,] -0.01647023314  0.7112123            0            0 0.0021333333
##  [16,] -0.00026800254  0.2458204            0            0 0.0014666667
##  [17,] -0.01152744851  0.3727665            0            0 0.0025333333
##  [18,] -0.00362124603  0.4932116            0            0 0.0025333333
##  [19,] -0.02240938075  0.8128801            0            0 0.0028000000
##  [20,]  0.01113662913  0.4821056            0            0 0.0016000000
##  [21,] -0.02836486280  0.8092080            0            0 0.0025333333
##  [22,] -0.10368334545  1.9930455            0            0 0.0058666667
##  [23,] -0.06378075968  1.4115292            0            0 0.0036000000
##  [24,] -0.17107318957  2.4115089            0            0 0.0077333333
##  [25,] -0.00986094243  0.5076520            0            0 0.0025333333
##  [26,] -0.14250513955  2.2457617            0            0 0.0062666667
##  [27,] -0.04747840190  1.2283511            0            0 0.0034666667
##  [28,]  0.00978619120  0.4129487            0            0 0.0010666667
##  [29,] -0.00351463226  0.3097896            0            0 0.0020000000
##  [30,] -0.00738293457  0.3842999            0            0 0.0012000000
##  [31,]  0.00486774189  0.3790955            0            0 0.0021333333
##  [32,]  0.03384653345  1.0808633            0            0 0.0028000000
##  [33,]  0.00654153167  0.3261972            0            0 0.0009333333
##  [34,]  0.00625761431  0.4010372            0            0 0.0017333333
##  [35,]  0.06025408379  1.5010229            0            0 0.0030666667
##  [36,]  0.01155648563  0.4596850            0            0 0.0020000000
##  [37,]  0.00177964328  0.4078900            0            0 0.0021333333
##  [38,] -0.00193298294  0.3107592            0            0 0.0009333333
##  [39,] -0.00573018811  0.4239365            0            0 0.0022666667
##  [40,]  0.00082028627  0.5543440            0            0 0.0022666667
##  [41,]  0.01204786388  0.4962271            0            0 0.0018666667
##  [42,] -0.02178836559  0.9648496            0            0 0.0029333333
##  [43,] -0.01684532461  0.9548668            0            0 0.0028000000
##  [44,] -0.00675761965  0.3395086            0            0 0.0028000000
##  [45,] -0.00613913712  0.4891281            0            0 0.0016000000
##  [46,] -0.00762962299  0.4389237            0            0 0.0016000000
##  [47,] -0.00663410436  0.3756887            0            0 0.0016000000
##  [48,]  0.00669872724  0.3964017            0            0 0.0017333333
##  [49,] -0.00034435030  0.1390577            0            0 0.0010666667
##  [50,]  0.01833404806  0.8046079            0            0 0.0029333333
##  [51,]  0.04158558967  1.4713239            0            0 0.0022666667
##  [52,]  0.01122000351  0.5265119            0            0 0.0014666667
##  [53,]  0.00617330700  0.3277792            0            0 0.0016000000
##  [54,] -0.00319363424  0.2111496            0            0 0.0014666667
##  [55,]  0.01188933563  0.5087609            0            0 0.0025333333
##  [56,]  0.00167179752  0.3188443            0            0 0.0018666667
##  [57,] -0.00953491420  0.3848878            0            0 0.0014666667
##  [58,]  0.00932028287  0.6625620            0            0 0.0021333333
##  [59,] -0.01507875291  0.5634662            0            0 0.0018666667
##  [60,] -0.00903040186  0.3625436            0            0 0.0021333333
##  [61,]  0.00009563128  0.1890493            0            0 0.0012000000
##  [62,] -0.00106960931  0.2316240            0            0 0.0014666667
##  [63,] -0.00010282415  0.4312542            0            0 0.0024000000
##  [64,] -0.00819730528  0.4278985            0            0 0.0017333333
##  [65,] -0.00271505765  0.5355645            0            0 0.0024000000
##  [66,]  0.00696295275  0.3505754            0            0 0.0025333333
##  [67,]  0.00886277979  0.6603891            0            0 0.0021333333
##  [68,]  0.01227569826  0.4671397            0            0 0.0021333333
##  [69,]  0.02017797812  0.6489430            0            0 0.0033333333
##  [70,]  0.03460275885  1.0235756            0            0 0.0025333333
##  [71,] -0.00531811777  0.3803452            0            0 0.0017333333
##  [72,] -0.02482249166  0.7397177            0            0 0.0025333333
##  [73,]  0.00136461852  0.3020055            0            0 0.0021333333
##  [74,] -0.01021762302  0.5697642            0            0 0.0025333333
##  [75,] -0.02306432047  0.7495557            0            0 0.0030666667
##  [76,] -0.02784400404  0.8632075            0            0 0.0030666667
##  [77,]  0.00790505529  0.5734143            0            0 0.0025333333
##  [78,] -0.00056101225  0.2453785            0            0 0.0021333333
##  [79,]  0.01096704290  0.4862492            0            0 0.0034666667
##  [80,]  0.00635613678  0.4227833            0            0 0.0012000000
##  [81,] -0.00040256004  0.2218477            0            0 0.0014666667
##  [82,]  0.00515834614  0.5372605            0            0 0.0020000000
##  [83,] -0.00007035086  0.4211832            0            0 0.0022666667
##  [84,]  0.00610592089  0.4112330            0            0 0.0026666667
##  [85,] -0.00296153989  0.4252058            0            0 0.0025333333
##  [86,]  0.00233603166  0.4785649            0            0 0.0018666667
##  [87,]  0.00950760536  0.4375328            0            0 0.0016000000
##  [88,] -0.00319741519  0.2553187            0            0 0.0012000000
##  [89,]  0.00279567956  0.2591166            0            0 0.0012000000
##  [90,] -0.01095113162  0.7119509            0            0 0.0018666667
##  [91,]  0.01754758253  0.7789657            0            0 0.0021333333
##  [92,] -0.00248372582  0.3042499            0            0 0.0018666667
##  [93,]  0.01056954134  0.6007759            0            0 0.0026666667
##  [94,]  0.01354972302  0.8107644            0            0 0.0020000000
##  [95,]  0.00610162425  0.2657757            0            0 0.0018666667
##  [96,]  0.00550132255  0.3047765            0            0 0.0008000000
##  [97,] -0.00328594477  0.3239822            0            0 0.0013333333
##  [98,]  0.00824007988  0.5886503            0            0 0.0025333333
##  [99,]  0.00377064657  0.3345422            0            0 0.0017333333
## [100,]  0.00268087629  0.4923001            0            0 0.0017333333
## [101,]  0.00361759755  0.2664276            0            0 0.0016000000
## [102,] -0.00125785535  0.1884398            0            0 0.0020000000
## [103,]  0.00277632737  0.4470401            0            0 0.0022666667
## [104,] -0.01170916706  0.7070968            0            0 0.0024000000
## [105,] -0.00950454374  0.5092766            0            0 0.0014666667
## [106,]  0.00991839707  0.4102977            0            0 0.0018666667
## [107,]  0.00760498155  0.7110995            0            0 0.0030666667
## [108,] -0.00717883556  0.5464990            0            0 0.0029333333
## [109,] -0.00824565656  0.6836161            0            0 0.0024000000
## [110,]  0.00751171461  0.3115849            0            0 0.0013333333
## [111,] -0.00055621021  0.2717771            0            0 0.0021333333
## [112,]  0.00084782172  0.3858959            0            0 0.0018666667
## [113,] -0.01374735320  0.5549056            0            0 0.0024000000
## [114,] -0.00030905418  0.4253479            0            0 0.0016000000
## [115,]  0.00044446584  0.3504592            0            0 0.0020000000
## [116,]  0.00892639020  0.3272707            0            0 0.0020000000
## [117,]  0.03063198668  0.8941522            0            0 0.0024000000
## [118,]  0.05853421103  1.3669008            0            0 0.0036000000
## [119,] -0.01799978263  0.5900237            0            0 0.0024000000
## [120,] -0.01381650457  0.6588336            0            0 0.0024000000
## [121,] -0.00801616383  0.4098493            0            0 0.0012000000
## [122,]  0.00231720140  0.2151038            0            0 0.0014666667
## [123,] -0.02841886236  0.8972683            0            0 0.0034666667
## [124,] -0.02196784884  0.8551798            0            0 0.0024000000
## [125,] -0.01820907688  0.6794758            0            0 0.0021333333
## [126,]  0.00191917596  0.4796340            0            0 0.0028000000
## [127,] -0.05492967323  1.1487590            0            0 0.0037333333
## [128,]  0.00075666617  0.4204341            0            0 0.0024000000
## [129,] -0.00192662113  0.3921484            0            0 0.0017333333
## [130,] -0.00894344962  0.5742215            0            0 0.0017333333
## [131,] -0.01273999395  0.5222802            0            0 0.0012000000
## [132,] -0.00101176390  0.6148757            0            0 0.0016000000
## [133,] -0.01271177986  0.7720026            0            0 0.0018666667
## [134,]  0.00198568213  0.8165910            0            0 0.0024000000
## [135,] -0.00614984337  0.5660170            0            0 0.0020000000
## [136,]  0.00977359504  0.6418780            0            0 0.0020000000
## [137,]  0.01685598053  0.5467248            0            0 0.0020000000
## [138,]  0.00478643032  0.6623801            0            0 0.0022666667
## [139,] -0.00148768633  0.1901172            0            0 0.0008000000
## [140,]  0.03684251473  1.0686060            0            0 0.0025333333
## [141,]  0.02237133934  0.7820852            0            0 0.0025333333
## [142,]  0.08334718299  2.2104629            0            0 0.0041333333
## [143,]  0.01670596339  0.7741902            0            0 0.0026666667
## [144,] -0.00646702378  0.3409709            0            0 0.0010666667
## [145,] -0.01787038539  0.6178606            0            0 0.0020000000
## [146,] -0.01210710502  0.5491503            0            0 0.0017333333
## [147,] -0.01766597014  0.7662021            0            0 0.0021333333
## [148,] -0.13610005814  2.1474472            0            0 0.0064000000
## [149,] -0.03446115081  0.9651905            0            0 0.0024000000
## [150,] -0.05639773821  1.2668824            0            0 0.0033333333
## [151,] -0.02135982801  0.6118892            0            0 0.0017333333
## [152,] -0.00711709735  0.2611309            0            0 0.0016000000
## [153,] -0.00271554287  0.1810662            0            0 0.0006666667
## [154,] -0.01551848705  0.6425994            0            0 0.0029333333
## [155,]  0.00231839316  0.4860063            0            0 0.0024000000
## [156,] -0.05187495157  1.1410308            0            0 0.0034666667
## [157,] -0.00092254702  0.1006495            0            0 0.0014666667
## [158,] -0.00349315349  0.2194937            0            0 0.0017333333
## [159,]  0.00638799743  0.3075779            0            0 0.0016000000
## [160,] -0.00319303611  0.5067947            0            0 0.0017333333
## [161,]  0.00187973045  0.2029551            0            0 0.0010666667
## [162,] -0.00083055311  0.4175919            0            0 0.0014666667
## [163,]  0.00722467650  0.4821929            0            0 0.0022666667
## [164,]  0.01699282389  0.7166051            0            0 0.0028000000
## [165,] -0.00236839311  0.4269201            0            0 0.0028000000
## [166,] -0.00704964679  0.2913808            0            0 0.0010666667
## [167,]  0.00231751600  0.2206739            0            0 0.0016000000
## [168,]  0.00163668294  0.2172364            0            0 0.0010666667
## [169,] -0.01664527911  0.6274140            0            0 0.0020000000
## [170,] -0.00899307613  0.3302830            0            0 0.0018666667
## [171,]  0.00280111222  0.2434321            0            0 0.0014666667
## [172,] -0.00081297778  0.6886768            0            0 0.0032000000
## [173,] -0.03431340239  1.1892105            0            0 0.0032000000
## [174,]  0.01020203271  1.3276948            0            0 0.0034666667
## [175,] -0.00306244740  1.1023503            0            0 0.0025333333
## [176,] -0.04166310150  1.5189742            0            0 0.0025333333
## [177,] -0.04865627529  1.3699596            0            0 0.0029333333
## [178,]  0.02272181788  1.1637376            0            0 0.0021333333
## [179,] -0.06124256069  1.8908378            0            0 0.0033333333
## [180,] -0.09466697984  2.0758565            0            0 0.0037333333
## [181,]  0.00520553456  1.0137846            0            0 0.0026666667
## [182,]  0.00038624600  0.2136883            0            0 0.0016000000
## [183,] -0.00775767903  1.4737020            0            0 0.0028000000
## [184,] -0.00737470658  0.6535756            0            0 0.0024000000
## [185,]  0.01983007156  0.5953587            0            0 0.0028000000
## [186,] -0.00808493610  0.5810289            0            0 0.0028000000
## [187,] -0.01034594011  0.6195578            0            0 0.0017333333
## [188,] -0.01815461465  0.9762409            0            0 0.0020000000
## [189,] -0.02674016229  1.1013299            0            0 0.0020000000
## [190,] -0.01556922389  1.2275186            0            0 0.0024000000
## [191,] -0.09440391118  2.2734959            0            0 0.0037333333
## [192,] -0.11541966605  2.8579405            0            0 0.0042666667
## [193,] -0.06505030255  1.8299151            0            0 0.0030666667
## [194,] -0.01912226918  0.9089227            0            0 0.0026666667
## [195,]  0.02722410942  2.6070112            0            0 0.0030666667
## [196,] -0.06220270214  1.9149677            0            0 0.0028000000
## [197,] -0.05534092807  2.1367940            0            0 0.0030666667
## [198,] -0.00754190509  0.4270707            0            0 0.0016000000
## [199,] -0.03940616146  1.7275463            0            0 0.0029333333
## [200,] -0.38716295306  5.5544717            0            0 0.0082666667
## [201,]  0.04122867996  1.2606507            0            0 0.0028000000
## [202,] -0.01060857689  0.4979391            0            0 0.0021333333
## [203,]  0.00162961556  0.4884729            0            0 0.0021333333
## [204,]  0.01036658606  1.1853954            0            0 0.0025333333
## [205,] -0.00338327579  1.1089221            0            0 0.0020000000
## [206,]  0.00435476674  1.3530703            0            0 0.0028000000
## [207,]  0.05356144388  2.6682333            0            0 0.0026666667
## [208,] -0.05672995038  2.5429039            0            0 0.0028000000
## [209,]  0.04243628741  1.7509027            0            0 0.0024000000
## [210,] -0.02046091494  0.7706384            0            0 0.0028000000
## [211,] -0.02405431603  1.7130633            0            0 0.0032000000
## [212,] -0.02422102191  0.9172878            0            0 0.0030666667
## [213,]  0.01101646059  1.4250734            0            0 0.0020000000
## [214,] -0.00063976827  0.3329916            0            0 0.0018666667
## [215,] -0.00043308647  0.7750810            0            0 0.0013333333
## [216,] -0.00164320983  0.4185761            0            0 0.0022666667
## [217,]  0.00744330952  0.4244182            0            0 0.0014666667
## [218,]  0.00155443218  0.4258177            0            0 0.0013333333
## [219,] -0.00215161030  0.6208105            0            0 0.0018666667
## [220,] -0.02766362758  0.9629175            0            0 0.0029333333
## [221,] -0.00556598363  0.7027008            0            0 0.0032000000
## [222,] -0.02051517955  1.1680229            0            0 0.0020000000
## [223,]  0.02136577403  2.0597233            0            0 0.0021333333
## [224,]  0.01420880434  0.9716482            0            0 0.0026666667
## [225,] -0.03947772384  1.0719718            0            0 0.0029333333
## [226,] -0.00990273047  0.7471768            0            0 0.0018666667
## [227,] -0.02368929216  2.1886330            0            0 0.0025333333
## [228,] -0.02507830105  0.8714135            0            0 0.0029333333
## [229,]  0.00191519565  0.4649630            0            0 0.0021333333
## [230,] -0.02356780432  0.7582051            0            0 0.0025333333
## [231,]  0.01969205293  1.0462867            0            0 0.0028000000
## [232,]  0.01062621297  0.6339565            0            0 0.0017333333
## [233,] -0.04767319049  1.1191493            0            0 0.0034666667
## [234,] -0.00624219615  0.5686079            0            0 0.0013333333
## [235,] -0.07297478542  1.7677986            0            0 0.0037333333
## [236,]  0.09308434432  3.1620765            0            0 0.0037333333
## [237,]  0.01367961526  0.8796297            0            0 0.0017333333
## [238,] -0.03100237972  1.1793412            0            0 0.0032000000
## [239,]  0.18407676227  4.1054085            0            0 0.0045333333
## [240,]  0.01925506764  1.4312522            0            0 0.0020000000
## [241,] -0.00332800596  0.7648545            0            0 0.0026666667
## [242,] -0.02034019953  1.0077404            0            0 0.0034666667
## [243,] -0.01151928020  0.9615075            0            0 0.0026666667
## [244,] -0.08346867671  1.8703253            0            0 0.0048000000
## [245,]  0.01068482020  0.8011117            0            0 0.0014666667
## [246,] -0.12704828267  2.0904095            0            0 0.0058666667
## [247,] -0.01582609695  0.8666925            0            0 0.0022666667
## [248,] -0.00664341987  0.5205422            0            0 0.0025333333
## [249,] -0.00768390617  0.3883377            0            0 0.0018666667
## [250,]  0.01747066456  0.6650394            0            0 0.0018666667
## [251,] -0.00964438519  0.5120514            0            0 0.0026666667
## [252,]  0.01272525738  1.7705447            0            0 0.0016000000
## [253,] -0.01949159208  0.7819220            0            0 0.0025333333
## [254,] -0.01046251581  1.0086585            0            0 0.0014666667
## [255,]  0.23065547591  6.2896192            0            0 0.0032000000
## [256,] -0.06016821592  2.2611802            0            0 0.0036000000
## [257,]  0.01317283119  0.7031169            0            0 0.0026666667
## [258,] -0.00460046878  0.4615307            0            0 0.0025333333
## [259,] -0.10953812678  4.3456181            0            0 0.0034666667
## [260,] -0.00779575824  0.5559891            0            0 0.0020000000
## [261,]  0.00383455567  1.2674115            0            0 0.0022666667
## [262,] -0.00273202064  0.2961800            0            0 0.0020000000
## [263,] -0.02486766639  1.5239684            0            0 0.0029333333
## [264,] -0.03221595933  1.4038361            0            0 0.0021333333
## [265,]  0.00584595257  0.6471392            0            0 0.0022666667
## [266,] -0.01288032877  0.5282785            0            0 0.0022666667
## [267,]  0.00958446685  0.4937792            0            0 0.0010666667
## [268,] -0.00051543682  0.8231363            0            0 0.0017333333
## [269,]  0.00296068457  0.4619030            0            0 0.0022666667
## [270,] -0.01110336188  1.0910313            0            0 0.0021333333
## [271,]  0.10569172239  3.2441184            0            0 0.0037333333
## [272,] -0.00500407114  0.7799976            0            0 0.0013333333
## [273,]  0.01348582480  1.0429696            0            0 0.0026666667
## [274,]  0.00214440035  1.0787527            0            0 0.0026666667
## [275,] -0.03961722952  1.8730041            0            0 0.0033333333
## [276,]  0.01139647608  0.8074230            0            0 0.0024000000
## [277,] -0.01103854317  0.6308667            0            0 0.0017333333
## [278,] -0.01205514139  0.6580385            0            0 0.0029333333
## [279,]  0.00257183737  1.2283511            0            0 0.0029333333
## [280,] -0.09736940558  2.4403745            0            0 0.0038666667
## [281,]  0.01163660297  0.4175403            0            0 0.0020000000
## [282,] -0.06649922248  1.5648285            0            0 0.0037333333
## [283,]  0.01098834837  0.8059703            0            0 0.0028000000
## [284,] -0.00052475651  0.6636149            0            0 0.0024000000
## [285,] -0.00718199968  0.9689494            0            0 0.0022666667
## [286,]  0.05676040446  2.5439325            0            0 0.0030666667
## [287,] -0.00077446957  1.2032435            0            0 0.0022666667
## [288,]  0.03798122587  1.6691974            0            0 0.0018666667
## [289,] -0.01888351315  0.7508118            0            0 0.0025333333
## [290,]  0.06999739199  1.8038485            0            0 0.0038666667
## [291,] -0.03747363234  1.5030872            0            0 0.0025333333
## [292,]  0.05936526073  1.9413802            0            0 0.0029333333
## [293,]  0.01886938863  1.4535352            0            0 0.0026666667
## [294,]  0.00174599762  0.4802649            0            0 0.0021333333
## [295,]  0.02146633259  1.4784986            0            0 0.0017333333
## [296,] -0.01003938128  0.8981341            0            0 0.0021333333
## [297,]  0.00246918798  0.4344574            0            0 0.0018666667
## [298,] -0.00214915651  0.4119082            0            0 0.0017333333
## [299,]  0.03585327099  1.4999970            0            0 0.0021333333
## [300,] -0.01204974994  0.7223358            0            0 0.0033333333
## [301,]  0.10001951896  3.3922676            0            0 0.0038666667
## [302,]  0.00795935423  0.6471900            0            0 0.0020000000
## [303,]  0.00756847533  1.0036607            0            0 0.0030666667
## [304,]  0.00462243457  0.8571579            0            0 0.0032000000
## [305,]  0.00618317451  0.7161878            0            0 0.0025333333
## [306,] -0.00308465839  0.8810382            0            0 0.0020000000
## [307,] -0.03995431964  2.9813258            0            0 0.0021333333
## [308,]  0.03570797249  1.2577194            0            0 0.0032000000
## [309,] -0.00590692519  0.4147279            0            0 0.0020000000
## [310,] -0.00587217822  0.2587566            0            0 0.0017333333
## [311,] -0.00419693225  0.8878103            0            0 0.0021333333
## [312,] -0.00358781651  0.7196047            0            0 0.0020000000
## [313,]  0.00009117405  0.2979145            0            0 0.0021333333
## [314,] -0.00794283133  0.6021553            0            0 0.0021333333
## [315,]  0.01845327463  0.8951495            0            0 0.0022666667
## [316,] -0.02301533698  0.9855940            0            0 0.0016000000
## [317,] -0.02322966334  0.9100376            0            0 0.0028000000
## [318,]  0.00878496706  0.8142470            0            0 0.0022666667
## [319,]  0.07016490094  2.9730100            0            0 0.0033333333
## [320,] -0.02754680319  0.9148212            0            0 0.0029333333
## [321,] -0.00405799753  0.6138842            0            0 0.0030666667
## [322,] -0.01381030436  0.9571022            0            0 0.0024000000
## [323,] -0.35308212308  9.4437342            0            0 0.0042666667
## [324,] -0.00695767588  0.6764322            0            0 0.0017333333
## [325,]  0.00374329079  0.3925999            0            0 0.0021333333
## [326,]  0.00018037313  0.6355554            0            0 0.0020000000
## [327,] -0.02371853914  0.9696055            0            0 0.0025333333
## [328,] -0.06570402085  1.7301026            0            0 0.0033333333
## [329,] -0.00375723919  0.2468508            0            0 0.0017333333
## [330,]  0.00094764973  0.2418656            0            0 0.0013333333
## [331,] -0.00119904665  0.6752365            0            0 0.0018666667
## [332,] -0.01255210481  0.8980907            0            0 0.0017333333
## [333,] -0.02762448463  0.9563711            0            0 0.0022666667
## [334,]  0.01351343700  1.7914729            0            0 0.0026666667
## [335,] -0.02080827693  0.9893235            0            0 0.0026666667
## [336,] -0.07854615426  2.0733521            0            0 0.0028000000
## [337,]  0.01313104320  0.7623468            0            0 0.0024000000
## [338,] -0.00973182484  0.4518883            0            0 0.0013333333
## [339,] -0.00206754392  1.0372468            0            0 0.0030666667
## [340,] -0.01999001313  0.9746472            0            0 0.0029333333
## [341,]  0.00071573863  0.3526947            0            0 0.0018666667
## [342,]  0.00260855587  0.4130259            0            0 0.0022666667
## [343,] -0.01078974107  0.9436231            0            0 0.0024000000
## [344,] -0.09997933831  2.4714020            0            0 0.0036000000
## [345,]  0.00581780556  0.5477047            0            0 0.0020000000
## [346,]  0.03464659912  0.9728667            0            0 0.0032000000
## [347,] -0.02746845970  0.8049390            0            0 0.0030666667
## [348,] -0.05196948929  1.6847849            0            0 0.0036000000
## [349,] -0.00770447005  0.9568034            0            0 0.0024000000
## [350,] -0.03548639673  1.3588053            0            0 0.0030666667
## [351,]  0.01746270507  2.0771395            0            0 0.0018666667
## [352,] -0.00976209087  0.5721004            0            0 0.0018666667
## [353,]  0.01307056949  1.1395545            0            0 0.0021333333
## [354,] -0.07073266794  2.1072539            0            0 0.0042666667
## [355,] -0.00327447574  0.5864519            0            0 0.0028000000
## [356,] -0.00159596872  0.6706603            0            0 0.0020000000
## [357,]  0.00504552329  0.8321624            0            0 0.0024000000
## [358,]  0.00854416068  0.4612553            0            0 0.0028000000
## [359,] -0.05274211642  2.0129653            0            0 0.0028000000
## [360,] -0.00377742376  0.6445630            0            0 0.0018666667
## [361,] -0.01158475698  0.6357252            0            0 0.0022666667
## [362,] -0.00008417772  0.4661515            0            0 0.0012000000
## [363,] -0.02499419521  0.9011908            0            0 0.0022666667
## [364,] -0.00443593796  0.3477457            0            0 0.0016000000
## [365,] -0.01049201805  0.5550310            0            0 0.0017333333
## [366,] -0.04586069885  1.4326092            0            0 0.0028000000
## [367,] -0.01109456341  0.5960701            0            0 0.0018666667
## [368,] -0.03920229652  1.3148421            0            0 0.0022666667
## [369,] -0.00129635499  0.5275291            0            0 0.0022666667
## [370,] -0.06081634329  1.8757021            0            0 0.0034666667
## [371,] -0.20830460427  5.8098532            0            0 0.0034666667
## [372,]  0.01657636722  0.9601838            0            0 0.0022666667
## [373,] -0.03899514392  1.7358351            0            0 0.0026666667
## [374,]  0.00365164172  0.6491723            0            0 0.0022666667
## [375,] -0.06984778067  2.0917160            0            0 0.0036000000
## [376,] -0.00452754283  0.7363184            0            0 0.0020000000
## [377,] -0.00537875812  0.3040461            0            0 0.0014666667
## [378,]  0.00344620705  0.2859794            0            0 0.0016000000
## [379,] -0.01373807322  0.5825774            0            0 0.0020000000
## [380,]  0.01395358919  1.1906888            0            0 0.0029333333
## [381,] -0.01840620307  1.0885049            0            0 0.0021333333
## [382,]  0.27126018601  8.1865463            0            0 0.0036000000
## [383,] -0.00070503954  1.2138388            0            0 0.0028000000
## [384,]  0.01373308238  1.0715478            0            0 0.0020000000
## [385,] -0.00267858361  0.7801755            0            0 0.0026666667
## [386,] -0.01676222307  0.6943828            0            0 0.0024000000
## [387,] -0.53665600360 12.2999086            0            0 0.0045333333
## [388,] -0.02291649049  0.9412741            0            0 0.0024000000
## [389,]  0.02127266261  0.9879992            0            0 0.0017333333
## [390,] -0.00346625202  0.6666998            0            0 0.0022666667
## [391,]  0.00132013458  1.6238668            0            0 0.0028000000
## [392,]  0.01196266853  0.8383648            0            0 0.0017333333
## [393,] -0.05426469693  1.3924142            0            0 0.0041333333
## [394,] -0.01865640157  0.5815308            0            0 0.0021333333
## [395,] -0.01084343685  0.6510940            0            0 0.0024000000
## [396,]  0.00323345240  0.5676648            0            0 0.0017333333
## [397,] -0.01327545731  0.7330531            0            0 0.0021333333
## [398,]  0.21992452966  9.1761354            0            0 0.0033333333
## [399,]  0.07591143163  5.3299688            0            0 0.0028000000
## [400,] -0.01186001476  0.9673010            0            0 0.0020000000
## [401,] -0.00361229405  0.6028003            0            0 0.0024000000
## [402,] -0.03123105574  1.0227672            0            0 0.0022666667
## [403,] -0.20809723888  6.9799452            0            0 0.0037333333
## [404,] -0.03460902861  1.3838552            0            0 0.0024000000
## [405,]  0.02214186052  1.1211032            0            0 0.0028000000
## [406,] -0.00603835458  0.5421071            0            0 0.0028000000
## [407,] -0.00324938157  0.3575859            0            0 0.0016000000
## [408,] -0.00626684148  0.8303860            0            0 0.0029333333
## [409,] -0.06835406643  1.5391718            0            0 0.0038666667
## [410,] -0.00690732290  0.3817615            0            0 0.0016000000
## [411,] -0.01980612285  0.9578651            0            0 0.0024000000
## [412,] -0.00673767907  0.5213549            0            0 0.0009333333
## [413,]  0.01489563367  0.8773850            0            0 0.0026666667
## [414,] -0.00180245024  0.4714237            0            0 0.0022666667
## [415,] -0.03341480599  1.0769096            0            0 0.0026666667
## [416,] -0.00326013021  0.6229288            0            0 0.0032000000
## [417,] -0.01017460564  0.9479096            0            0 0.0033333333
## [418,] -0.01942328620  0.7012808            0            0 0.0022666667
## [419,]  0.01150827917  1.5981117            0            0 0.0021333333
## [420,] -0.00872144185  0.7176251            0            0 0.0020000000
## [421,] -0.02406431606  0.9412862            0            0 0.0024000000
## [422,] -0.01574707763  0.9365421            0            0 0.0026666667
## [423,] -0.01555830999  1.7705418            0            0 0.0025333333
## [424,] -0.00005970490  0.6676579            0            0 0.0024000000
## [425,] -0.00230019706  0.4073801            0            0 0.0014666667
## [426,] -0.00147609517  0.2752254            0            0 0.0022666667
## [427,] -0.01950261775  0.8638737            0            0 0.0026666667
## [428,] -0.01425060103  0.6446510            0            0 0.0025333333
## [429,] -0.00705720395  0.5807070            0            0 0.0026666667
## [430,] -0.02623462985  0.8793396            0            0 0.0028000000
## [431,] -0.07339784456  1.9810074            0            0 0.0033333333
## [432,]  0.03737303670  1.3511437            0            0 0.0029333333
## [433,] -0.02254038964  0.9771077            0            0 0.0032000000
## [434,] -0.10306606179  2.8014589            0            0 0.0049333333
## [435,] -0.00839785496  1.3698286            0            0 0.0025333333
## [436,] -0.04891872454  1.3471521            0            0 0.0030666667
## [437,] -0.01298447161  0.8422718            0            0 0.0030666667
## [438,] -0.00197269586  0.3240561            0            0 0.0012000000
## [439,] -0.04287659738  1.5395138            0            0 0.0033333333
## [440,] -0.05220001233  1.6482485            0            0 0.0033333333
## [441,] -0.00850247362  0.5605286            0            0 0.0021333333
## [442,] -0.00638144708  0.5194437            0            0 0.0025333333
## [443,] -0.12408356856  2.7675065            0            0 0.0042666667
## [444,] -0.00905194322  0.7470167            0            0 0.0029333333
## [445,] -0.00812968857  0.7944392            0            0 0.0022666667
## [446,]  0.00001768554  1.0378946            0            0 0.0028000000
## [447,] -0.02218304663  0.9107648            0            0 0.0028000000
## [448,]  0.00092772833  0.6332493            0            0 0.0021333333
## [449,] -0.02568873915  0.9461023            0            0 0.0030666667
## [450,]  0.01392776984  1.4068940            0            0 0.0020000000
## [451,] -0.06793815607  2.2470119            0            0 0.0032000000
## [452,]  0.00436084773  0.7823230            0            0 0.0032000000
## [453,]  0.00366257960  0.4327661            0            0 0.0016000000
## [454,]  0.00268852001  0.4048095            0            0 0.0021333333
## [455,] -0.00480479394  1.3534395            0            0 0.0029333333
## [456,] -0.02585450874  1.0825693            0            0 0.0020000000
## [457,] -0.00902950440  0.4714410            0            0 0.0021333333
## [458,] -0.01651540909  0.6028711            0            0 0.0020000000
## [459,] -0.00387793363  0.3980457            0            0 0.0021333333
## [460,] -0.00002073313  0.3378363            0            0 0.0021333333
## [461,]  0.00499728073  1.1117109            0            0 0.0032000000
## [462,]  0.00794450185  1.0856762            0            0 0.0028000000
## [463,] -0.02138763373  1.2370726            0            0 0.0032000000
## [464,]  0.01974010483  0.9974529            0            0 0.0026666667
## [465,]  0.00483055702  0.4848990            0            0 0.0021333333
## [466,] -0.00646321734  0.5308886            0            0 0.0021333333
## [467,] -0.10459797360  4.1466932            0            0 0.0033333333
## [468,]  0.01493883877  1.0473745            0            0 0.0021333333
## [469,]  0.02249751131  0.9571220            0            0 0.0025333333
## [470,] -0.00837668430  0.4257489            0            0 0.0017333333
## [471,] -0.03773395461  1.0908663            0            0 0.0033333333
## [472,]  0.01837954520  0.8211637            0            0 0.0022666667
## [473,] -0.04562452911  1.0408843            0            0 0.0030666667
## [474,] -0.01371617775  0.8715611            0            0 0.0022666667
## [475,] -0.06018091995  1.5702828            0            0 0.0037333333

4.2.3 Predict fitted values for each individual

pred.npb <- predict(fit.npb)
fittedvals <- pred.npb$fitted.vals

4.2.4 Plot predicted outcomes against “measured” outcomes

plot(fittedvals, Y)
abline(a = 0, b = 1, col = "red")